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Abstract 


Many  studies  of  aircraft  high  angle  of  attack  dynamics  have  argued  the  need 
for  rotary  aerodynamic  data  to  be  included  in  a  model  in  order  to  effectively  model 
aircraft  spin  behavior.  The  purpose  of  this  research  paper  was  to  use  bifurcation 
analysis  to  investigate  the  effectiveness  of  rotary  balance  data  in  the  prediction  of 
aircraft  spin  behavior  as  both  a  stand  alone  representation  of  a  model’s 
aerodynamic  data  and  in  a  conventional  hybrid  model.  Equilibrium  solutions  from 
both  models  were  compared  to  previous  studies  which  utilized  conventional  static 
and  forced  oscillation  aerodynamic  data  and  to  flight  test  results  to  analyze  the 
effectiveness  of  rotary  aerodynamic  data  for  the  prediction  of  spin  behavior. 

Using  the  foundation  of  a  previously  developed  model  of  the  F-15B,  the  rotary 
balance  aerodynamic  force  and  moment  data  were  implemented  as  a  function  of 
angle  of  attack,  angle  of  sideslip,  non-dimensional  rotation  steady  state  rate  and 
the  control  surface  deflections  (§a,5^j,5e,5,).  Bifurcation  diagrams  were  developed 
as  a  function  of  alpha  versus  and  Sg  to  show  highlights  of  equilibrium  and 
dynamic  behavior  of  the  aircraft.  For  selected  configurations,  the  resulting  aircraft 
state  variables  showed  the  rotary  balance  data  model  having  close  correlation  to 
experimental  flight  test  data.  Corhparison  of  these  selected  configurations  with  the 
hybrid  and  conventional  static  and  forced  oscillation  models,  showed  comparable 
results.  However,  the  models  bifurcation  diagrams  were  very  different.  Problems 
were  identified  with  static  contributions  of  the  rotary  balance  data  indicating  a 


Xli 


possible  cause.  Despite  the  static  contribution  problem,  the  overall  results 
indicated  the  rotary  balance  data  model  did  provide  a  reasonable  representation 
of  the  spin  dynamics  of  the  aircraft.  The  development  of  the  hybrid  model 
displayed  the  difficulties  in  blending  of  the  aerodynamic  coefficient  data  in  the 
presence  of  deficient  experimental  data,  inaccurate  modelling  of  aerodynamic 
coefficients  or  possible  differences  in  the  databases  such  as  Reynolds  number 
effects.  Recommendations  on  continued  investigation  of  the  effects  of  the  static 
contributions  of  the  rotary  balance  data  were  made. 


INVESTIGATION  OF  THE  INFLUENCE  OF 
ROTARY  AERODYNAMICS  ON  THE  STUDY  OF 
HIGH  ANGLE  OF  ATTACK  DYNAMICS  OF  THE 
F-15B  USING  BIFURCATION  ANALYSIS 


I.  Introduction 


The  flat  spin  is  the  most  dangerous  spinning  motion  exhibited  by  an  airplane 
where  the  aircraft’s  angle  of  attack  (AOA)  approaches  90°  with  a  high  rate  of 
rotation.  Conventional  aircraft  control  surfaces  become  increasingly  ineffective  as 
the  AOA  approaches  90°  and  as  the  spin  develops  it  may  become  impossible  for 
the  pilot  to  recover.  Today’s  fighter  aircraft  are  also  heavier  than  earlier  designs, 
with  weight  more  concentrated  along  the  fuselage.  This  design  evolution  was  a 
factor  of  the  desirability  of  a  more  maneuverable  aircraft.  With  this  distribution  of 
mass,  the  yaw  moment  of  inertia  of  the  aircraft  has  increased  as  much  as  20  times 
compared  to  early  fighters  but  the  effective  control  surfaces  have  remained 
practically  the  same  (11:1).  This  combination  of  high  yaw  moment  of  inertia  and 
relatively  ineffective  control  surfaces  at  high  angles  of  attack  warrants  analysis  of 
spin  characteristics  of  aircraft  because  recovery  is  very  difficult  if  not  fatal. 

Aircraft  high  angle-of-attack  research  has  been  increasing  over  the  past  twenty 
years.  The  most  recent  full  scale  research  effort  is  the  X-31  where  one  of  its 
areas  of  investigation  is  the  use  of  thrust  vectoring  paddles  to  improve  the  agility 
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and  handling  qualities  of  an  aircraft  while  flying  at  high  angles  of  attack.  In 
February  1991,  testing  began  at  AOA  up  to  70°.  The  ability  of  an  aircraft  to 
achieve  high  angle  of  attack  will  considerably  reduce  the  time  needed  to  maneuver 
an  aircraft  to  obtain  a  good  firing  position  in  close-in  combat  (14:38).  Research 
is  ongoing  to  improve  designs  of  future  aircraft  and  their  control  systems  to  have 
stability  and  agility  during  high  AOA  maneuvers.  For  developed  aircraft,  it  is 
desired  to  design  either  effective  pilot  procedures  or  aircraft  control  systems  that 
can  predict  and  recover  an  aircraft  from  potentially  fatal  attitudes  accompanied  with 
high  angle  of  attack  maneuvers.  Loss  of  control  of  the  aircraft  can  occur  through 
non-linear  behavior  such  as  stalls,  departures,  wing  rock,  nose  slice,  spin  entry, 
and  full  spins.  Considering  the  costs  involved  in  full  scale  testing,  it  is  desirable 
for  the  development  of  a  methodology  that  could  investigate  the  aircraft  dynamics 
associated  with  high  angle  of  attack  flight,  specifically  aircraft  spin  behavior,  with 
the  use  of  scale  model  aerodynamic  data.  Investigations  of  this  form  have  used 
combinations  of  static  wind  tunnel,  forced  oscillation  and  rotary  aerodynamic  data. 
It  is  the  objective  of  this  paper  to  investigate  the  effectiveness  of  rotary 
aerodynamic  data  in  modeling  high  angle  of  attack  dynamics  in  the  regime  of 
a=30°  to  90°  using  the  methodology  of  bifurcation  analysis. 

Previous  Studies 


The  mathematical  modeling  of  an  aircraft’s  motion  during  high  angle  of  attack 
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maneuvering  has  shown  to  be  an  effective  tool  for  the  flight  dynamicist  however 
it  is  highly  dependent  on  the  results  of  wind  tunnel  tests  for  the  required  static  and 
dynamic  aerodynamic  data  (27:244).  In  1972,  Adams  (1 )  used  a  numerical  routine 
to  predict  spin  modes  for  various  aircraft  by  searching  for  steady  spin  states.  His 
results  did  not  compare  well  with  flight  test  data.  Adams  attributed  poor  results  to 
the  deficiencies  of  the  aerodynamic  data  his  model  was  based  on.  The 
discrepancies  he  found  may  have  been  corrected  with  the  inclusion  of  rotary 
balance  data. 

As  technology  advanced  with  digital  computers  in  the  1960s  many  flight 
dynamic  models  of  aircraft  spinning  were  developed.  In  most  cases  these  models 
involved  the  combination  of  aerodynamic  data  from  static  tests  with  small 
disturbance  oscillatory  data.  In  these  models  the  aerodynamic  oata  were 
sometimes  inaccurate  and  spin  modes  could  not  be  predicted  (22;'*44).  Early 
studies  of  spins  involved  numerical  simulations  using  static  wind  tunnel  data  to 
predict  spin  entry  and  possible  recovery  techniques.  As  indicated  above,  their 
results  did  not  compare  well  with  flight  test  data  nor  spin  tunnel  results. 
Chambers,  Bauman  and  Anglin  determined  that  rotary  balance  data  was  necessary 
to  correctly  model  the  aerodynamics  during  a  spin  (11).  in  many  investigations 
utilizing  rotary  balance  data,  this  conclusion  has  been  substantiated. 

During  the  past  forty  years  many  research  efforts  have  been  made  on  the  study 
of  aircraft  spin  behavior  utilizing  rotary  aerodynamic  data.  Some  encountered 
problems  obtaining  acceptable  results  utilizing  the  rotary  balance  data.  However 
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most  have  found  significance  of  including  rotary  aerodynamic  data  in  their  models 
for  prediction  of  aircraft  spin  modes.  In  1954,  Scher  (32)  applied  equations  of 
motion  and  spin-geometry  relationships  with  rotary  balance  aerodynamic  data  to 
calculate  step  by  step  the  details  of  a  spin  recovery  motion  for  an  unswept-wing 
fighter-airplane  configuration.  At  first,  his  results  did  not  agree  with  the 
assumptions  of  steady  spin.  However,  as  also  determined  by  Stone,  Burk  and 
Senger  (37),  the  rotary  balance  data  had  to  be  modified  due  to  inconsistencies  in 
the  mounting  of  the  model.  The  values  he  obtained  were  inconsistent  with  the 
assumption  that  the  aerodynamic  and  inertial  forces  and  moments  should  balance 
during  steady  spin.  As  a  result,  the  spin  behavicr  he  obtained  were  not  validated 
by  ttie  test  aircraft  of  the  investigation. 

With  the  improvement  of  the  technique  of  rotar^  balance  wind  tunnel  testing,  the 
integrity  of  the  data  has  improved  from  the  tests  of  over  forty  years  ago.  A  most 
recent  analysis  was  performed  in  1989  by  Martin  and  Hill  (22)  where  rotary 
balance  data  was  used  in  a  model  of  a  basic  training  aircraft  using  a  six-degree 
of  freedom  flight  dynamics  model  of  aircraft  spin.  The  model  was  used  to  predict 
equilibrium  spin  conditions  and  spin  recovery  techni...v'r,5.  Their  results  were 
promising  when  compared  to  scale  model  wind  tunnel  tests.  In  1981 ,  Tischler  and 
Barlow  (38)  were  able  to  accurately  determine  spin  modes  of  a  low-wing  general 
aviation  aircraft  using  rotary  balance  data.  Studies  by  Birhie  (9)  on  a  fighter 
aircraft  design  in  1980,  had  demonstrated  that  the  spin  computed  with  static 
aerodynamic  data  did  not  match  the  flight  motion  whereas  the  spin  computed  with 


rotational  data  duplicated  the  developed  spin.  This  result  again  emphasized  the 
significance  of  rotary  balance  data  representing  the  dynamics  during  spin  motion, 
Ogburn,  Nguyen  and  Hoffler  (27)  have  dr.iu'  • -.cent  research  showing  that  the 
“  addition  of  dynamic  terms,  including  ro ae  'Cv  amics  can  significantly  affect 
the  simulated  flight  motions.  Their  efforts  al;  o  showed  a  substaiitial  influence  of 
the  rotary  aerodynamic  data  on  the  mpi'  ,-  and  cont-'ols  for  the  aircraft 
configurations  they  tested. 

High  angle  of  attack  maneuvers  are  very  nonlinear  and  require  complex 
analysis  tools  to  study  the  behavior,  Mehra  and  Carroll  (26)  performed  fairly 
extensive  analysis  of  the  F-4  Phantom  fighter  in  1979  demonstrating  the  use  of 
bifurcation  and  continuation  methodology  as  an  analysis  tool  for  the  study  of 
aircraft  high  AOA  behavior.  Bifurcation  theory  was  not  a  new  concept  however  it 
had  not  been  previously  applieu  lO  aircraft  dynamics.  The  tool  enables  an  analyst 
to  use  more  complex  (higher  order,  nonlinear)  aircraft  models  thereby  enabling 
analysis  of  more  demanding  flight  conditions  such  as  spins,  stalls  and  wing  rock. 
Bifurcation  methodology  provides  insight  into  the  solution  of  nonlinear  equations 
through  development  of  a  mapping  of  equilibrium  and  periodic  solutions  of  an 
aircraft’s  equations  of  motion.  This  results  in  a  global  view  of  the  nonlinear 
behavior  of  the  aircraft.  Further  discussion  of  bifurcation  theory  is  presented  in 
Chapter  ill. 

Many  recent  research  efforts,  have  used  bifurcation  analysis  for  the  study  of 
aircraft  dynamics,  Hawkins  (17),  Jahnke  (18,19,20),  Zagaynov  and  Goman  (39), 
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and  Guicheteau  (16}  used  bifurcation  theory  to  analyze  nonlinear  behavior 
associated  \with  high  AOA  maneuvers  of  various  aircraft  configurations  including 
the  F-14  and  F-15  fighters.  Barth  (G)  and  Planeaux  and  Barth  (28)  did  an 
investigation  of  high  AOA  behavior  of  the  F-15  using  bifurcation  analysis  ar-d  were 
able  to  identify  periodic  solutions  indicative  of  aircraft  wing  rock  behavior. 
Baumann  (7),  Beck  (8)  and  Planeaux,  Beck  and  Baumann  (29)  continued  the  F-15 
research  with  an  expanded  aerodynamic  database  allowing  analysis  up  to  90° 
AOA.  With  the  new  model  they  analyzed  the  effectiveness  of  control 
augmentation.  The  most  recent  effort  with  the  F-15  model,  performed  by 
McDonnell  (24)  and  Planeaux  and  McDonnell  (30)  investigated  the  effectiveness 
of  thrust  vectoring  for  spin  recovery.  Many  options  for  spin  recovery  using  thrust 
vectoring  were  identified. 

Of  the  efforts  mentioned  utilizing  bifurcation  methodology,  only  three  included 
the  contributions  of  rotary  aerodynamic  data.  Hawkins  (17),  Jahnke  (28)  and 
Mehra  and  Carroll  (26,1 0)  developed  hybrid  models  combining  the  data  from  static, 
forced  c<sci!lation  and  rotary  balance  testing  in  attempt,  to  more  accurately  model 
the  aircraft  dynamics.  Mehra  and  Carroll’s  analysis  included  an  investigation  of 
different  approaches  to  integrating  rotary  aerodynamic  data  into  a  model.  Their 
effort  did  provide  many  methods  for  consideration  some  of  which  were 
incorporated  by  Hawkins  and  Jahnke  in  their  aircraft  models. 
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Overview 


There  is  still  a  debate  over  what  effective  benefit  rotary  balance  data  provides 
in  a  model.  This  paper  will  continue  the  development  of  the  F-15B  model  by 
Barth  (6),  Beck  (8),  Baumann  (7)  and  McDonnell  (24).  The  inclusion  of  rotary 
balance  data  will  be  analyzed  in  two  developments.  The  rotary  balance  data  will 
first  be  used  in  the  model  as  the  principal  representation  of  the  static,  oscillatory 
and  rotary  aerodynamic  force  and  moment  coefficients.  The  results  of  bifurcation 
analysis  and  continuation  theory  will  be  used  to  compare  the  static  and  forced 
oscillation  aerodynamic  based  model  and  a  stand  alone  rotary  aerodynamic  based 
model.  The  comparison  will  assist  in  identification  of  the  strength  and  weakness 
of  rotary  balance  data  as  well  as  possible  deficiencies  in  the  modeling  of  the 
coefficient  data.  The  second  phase  of  the  investigation  will  be  the  development 
and  analysis  of  a  hybrid  model  utilizing  the  rotary  aerodynamic  data  through 
blending  with  the  static  and  forced  oscillation  aerodynamic  data  based  model.  In 
addition,  based  on  problems  encountered  with  the  rotary  balance  database,  an 
analysis  of  static  aspects  of  the  rotary  balance  database  and  perturbation  of  two 
Oi  the  rotary  balance  coefficients  will  be  made. 

Chapter  II  will  discuss  spin  tunnel  testing,  the  rc^r^  balance  data  used  in  this 
analysis  and  the  data  processing  required.  Chapter  III  will  present  a  brief 
introduction  to  bifurcation  theory  and  the  continuation  methodology.  Chapter  IV 
will  describe  the  F-15  model  and  the  modifications  required  with  implementation 
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of  the  rotary  balance  data.  Chapter  V  presents  the  findings  of  the  analysis.  The 
conclusions  will  be  drawn  as  well  as  recommendalions  for  future  research  in 
Chapter  VI. 
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11.  Rotary  Balance  Data 

It  should  be  rather  intuitive  that  the  best  method  to  investigate  a  particular 
motion  of  an  aircraft  would  be  to  use  a  technique  that  most  accurately  represents 
that  motion.  The  standard  wind  tunnel  aerodynamic  coefficient  data  consists  of 
two  types,  static  and  forced  oscillation.  Both  sets  of  data  have  been  used  as  the 
foundation  for  dynamicist  to  study  the  dynamics  of  aircraft.  For  the  study  of 
aircraft  in  spinning  or  coning  motion,  the  most  representative  modeling  would  be 
to  rotate  a  model  with  its  spin  axis  parallel  to  the  free  stream  velocity  of  the 
relative  wind.  This  reasoning  has  been  identified  for  the  study  of  aircraft  research 
since  the  pioneer  days  of  aviation.  In  response,  rotary  balance  techniques  have 
been  developed  providing  information  on  the  effects  of  rotational  rates  on  the 
overall  aerodynamic  forces  and  moments  acting  on  the  aircraft. 

The  fundamental  problem  of  the  dynamicist  is  to  properly  represent  the 
equations  of  motion  with  the  aerodynamic  reactions  and  motion  variables  (control 
parameters  and  states)  of  interest.  When  analyzing  low  angle  of  attack  dynamics, 
adequate  dynamic  predictions  can  be  made  using  a  linearized  small  perturbation 
model  with  constant  aerodynamic  derivatives.  However,  when  motions  involve 
large  variation  in  any  of  the  state  variables,  change  in  some  of  the  aerodynamic 
derivatives  may  vary  appreciably  with  factors  of  the  flight  conditions  thereby  the 
analysis  must  involve  non-linearities  to  properly  represent  the  dynamics.  It  is 
argued  that  use  of  the  combination  of  aerodynamic  data  of  oscillatory  and  rotary 
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motion  can  produce  the  data  required  for  calculation  of  flight  motions  including 
fairly  large  excursion  motions  such  as  high  angle  of  attack  departure  and  spin 
entry.  The  difficulties  involved  in  performing  such  tests  and  receiving  substantial 
data  concern  the  correlation  between  wind-tunnel  dynamic  test  data  which  are 
usually  available  from  separate  oscillatory  and  rotary  tests,  where  the  model  or 
its  motion  did  not  exactly  simulate  the  relevant  flight  motion  of  interest  (2:209). 

For  rotary  tests,  the  axis  of  rotation  is  parallel  to  the  wind  axis.  The 
aerodynamic  loads  are  measured  by  a  strain  gauge  aligned  along  the  model  body 
axis.  The  model  has  a  fixed  orientation  of  angle  of  attack,  angle  of  sideslip  and 
selected  control  deflections  to  the  wind  direction  as  the  model  is  rotated.  The 
model  is  seeing  a  constant  rate  of  motion  resulting  in  a  consistent  collection  of 
data.  In  oscillatory  tests  the  motion  is  not  of  constant  rate  thereby  the  flow  is 
constantly  changing  and  under  certain  conditions  the  measurements  of  the 
derivative  data  may  be  a  function  of  the  history  of  the  flow  (2:209).  At  low  angles 
of  attack  and  sideslip  the  flow  is  attached  to  virtually  all  of  the  surface  of  the  model 
and  fairly  good  agreement  of  the  derivative  data  is  found  between  the  rotary  and 
oscillatory  test  results.  However,  at  high  angle  of  attack,  there  is  increased  flow 
separation  and  lag  effects  that  may  introduce  significant  differences  between  the 
two  testing  results.  As  indicated  in  the  AGARD  Advisory  Report  (2:210),  tests 
comparing  rotary  and  oscillatory  data  have  attributed  this  lag  effect  for  their 
apparent  differences.  These  effects  may  influential  in  the  development  of  the 
hybrid  model.  These  effects  will  be  considered  when  comparing  the  stand  alone 
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rotary  balance  model  with  the  static  and  forced  oscillation  data  based  model. 

The  above  discussion  brings  up  the  issue  of  the  integrity  of  the  F-15  1988 
Aerobase  which  is  the  source  of  aircraft  aerodynamic  force  and  moment 
coefficients  for  the  McDonnell  (24)  and  Baumann  (7)  theses.  Part  of  this 
investigation  involves  both  the  blending  and  comparing  of  the  dynamics  resulting 
from  both  the  rotary  balance  and  the  1988  Aerobase  which  represents  the 
oscillatory  testing  data.  The  1988  F-15  Aerobase  was  obtained  from  the  F-15 
System  Program  Office  (SPO)  at  Wright-Patterson  AFB,  Ohio  and  has  been  used 
for  simulation  testing  of  the  aircraft.  This  database  is  not  well  documented  on  the 
exact  source  for  its  content.  However,  the  McDonnell  Douglas  Documentation  (23) 
indicates  that  the  stability  and  control  derivatives  were  derived  by  analysis  of  flight 
test  data  and  where  appropriate,  superseded  relative  v/ind  tunnel  based  estimates. 
Additionally,  in  order  to  provide  the  maximum  amount  of  range  on  flight  conditions, 
data  from  flight  test  results  were  utilized  to  adjust  the  associated  wind  tunnel 
coefficient  data  (23).  Neither  the  method  utilized  for  blending  the  flight  test  data 
nor  a  discussion  of  the  methods  of  data  processing  used  on  the  raw  wind  tunnel 
data  were  identified  in  the  documentation.  It  should  be  noted  that  Inconsistencies 
may  be  present  for  this  investigation  In  the  effectiveness  of  rotary  balance  data, 
when  the  baseline  model  to  be  compared  and  blended  consisted  of  actual  flight 
test  data.  The  1988  F-15  Aerobase  is  not  a  true  representation  of  the  limitations 
of  conventional  wind  tunnel  testing  results.  This  factor  v/ill  be  considered  when  the 
results  are  analyzed  in  Chapter  V. 
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Description  of  the  Rotary  Aerodynamic  Data 


The  rotary  balance  technique  measures  the  6  aerodynamic  force  and  moment 
coefficients  C^,  Cy,  C^,  Q,  C^,  and  as  a  function  of  the  3  state  variables:  spin 
rate,  angle  of  attack  and  angle  of  sideslip;  and  the  control  deflections:  6e,  5„ 

The  rotary  balance  data  obtained  for  this  analysis  is  documented  in  NASA 
Contractor  Report  3478  (5)  and  3479  (4).  The  database  was  obtained  from  the 
NASA  Langley  Spin  Tunnel  in  the  form  of  data  files  on  floppy  disk.  The  data  files 
were  arranged  by  aircraft  configuration  representative  of  the  data  presented  in 
reference  (5)  and  consisted  of  the  six  aerodynamic  force  and  moment  coefficients 
tabulated  with  their  associated  state  variables  and  control  deflections.  The  spin 
tunnel  data  consisted  of  configurations  ranging  from  the  build-up  of  individual 
airplane  components  (body .wings, tail),  the  basic  airplane  configured  with  various 
control  deflections  and  the  airplane  configured  with  conformal  fuel  tanks.  Since 
this  analysis  was  restricted  to  a  basic  F-1 5B  with  no  stores  carriage  nor  conformal 
fuel  tanks,  only  16  configurations  were  used  from  the  available  database,  table 
I  shows  the  16  configurations  that  were  used  to  represent  the  rotary  balance  data 
in  this  analysis. 

The  Langley  Spin  Tunnel  tests  were  conducted  at  a  free  stream  velocity  of  25 
ft/sec,  which  correspond  to  a  Reynolds  number  of  approximately  21 1 ,000  based 
on  the  model  wing  cord  of  1.33  ft.  Studies  by  G.N.  Malcolm  at  NASA-Ames 
(2:1 17)  on  the  sensitivity  of  the  rotational  flow  fields  to  Reynolds  number  variation. 
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Table  I  Rotary  Balanc'=‘  Database  Testing  Configurations 


Configuration  6^*  6^*  63*  5r* 


1 

0 

0 

0 

0 

0 

2 

0 

0 

6 

0 

0 

3 

0 

0 

11 

0 

0 

4 

0 

0 

6 

20 

0 

5 

0 

0 

-6 

-20 

30 

6 

0 

0 

6 

20 

-15 

7 

0 

0 

-11 

-20 

30 

8 

0 

-25 

0 

0 

0 

9 

-10 

0 

0 

0 

0 

10 

10 

0 

0 

0 

0 

11 

10 

0 

6 

0 

0 

12 

10 

0 

11 

0 

0 

13 

10 

0 

-6 

-20 

30 

14 

10 

0 

6 

20 

-15 

15 

10 

0 

-11 

-20 

30 

16 

10 

-25 

0 

0 

0 

Each  configuration  is 

tested  for 

all 

combinations  of 

a  =  8,10, 

12,14,16,18,20 

,25,30,35,40,45 

,50,55,60 

,70,80,90 

degrees 

and  |a.b/2V„| 

II 

0 

0 

0 

M 

/0.2, 

.0.3, 0.4 

(some  for 

0.5, 0.6, 

0.7,0. 9)  for 

clockwise 

and 

counter 

clockwise 

rotations . 

*  in  degrees 


have  shown  that  variation  in  Reynolds  number  does  have  a  significant  effect  on 
the  behavior  of  the  aerodynamic  derivative  coefficients.  Aircraft  with  slender  noses 
or  slender  forebodies  can  experience  large  side  forces  due  to  asymmetric 
separation  and  vortices  on  the  leeward  side  of  thi-:  oody.  It  has  been  found  that 
the  behavior  of  these  vortices  and  flow  separation  change  as  Reynolds  number 
varies.  The  high  Reynolds  number  effects  from  the  experimental  flight  test  data 
inherent  in  the  1988  Aerobase  may  introduce  noticeable  differences  in  the 
aerodynamic  coefficients  when  compared  to  the  rotary  balance  data.  Since  high 
Reynolds  number  rotary  balance  data  was  unavailable,  this  investigation  was 
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constrained  to  blending  of  the  dissimilar  Reynolds  number  data.  The  Reynolds 
number  effect  may  account  for  gross  differences  when  comparing  the  coefficient 
data  from  each  model.  For  all  tests,  the  spin  axis  passed  through  the  full  scale 
airplane  nominal  center  of  gravity  (CG)  location  for  AOA  above  30°.  Testing  below 
30°  AOA  required  shifting  of  the  spin  axis  to  allow  the  model  a  full  range  of  motion. 
Use  of  this  data  would  require  adjustment  of  the  coefficient  data  to  the  nominal  eg 
location.  For  this  analysis  only  a  >  30°  was  used.  The  McDonnell  (24)  analysis 
was  based  upon  an  aircraft  at  0.6  Mach  and  20,000  ft  altitude. 

The  data  was  presented  for  an  angle-of-attack  range  of  8°  to  90°,  and  clockwise 
and  counter-clockwise  rotations  covering  range  from  0  to  0.4.  Some 

configurations  were  presented  with  an  extended  range  from  0  to  0.9.  The 

data  used  for  this  investigation  is  not  the  specific  data  published  in  reference  (5) 
but  a  second  run  collection  of  the  same  configurations  performed  in  June-August 
1981.  The  integrity  of  the  axial  and  side  force  coefficient  data  may  be 
questionable  (40).  Analysis  of  the  data  files  identified  this  problem  and  will  be 
discussed  in  the  next  section.  A  complete  description  of  the  testing  apparatus  and 
model  are  in  reference  (5). 

As  shown  by  Table  I,  the  depth  of  the  rotary  balance  database  is  limited  by  the 
few  airplane  configurations  tested  for  the  basic  F-15.  This  investigation  will  limit 
itself  to  analysis  of  only  configurations  that  fit  within  these  restrictions.  Table  II 
shows  the  F-15B  operational  control  deflection  limitations  as  compared  to  the 
rotary  balance  database  limitations.  Unlike  the  static  and  oscillatory  test  data,  the 
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rotary  balance  data  will  not  be  adjusted  to  the  limitations  of  the  operational  flight 
configurations  and  will  thereby  provide  a  better  indication  of  the  capabilities  of  this 
data  as  a  stand  alone  representation  of  the  aircraft  dynamics. 


Table  II  Control  Deflection  Limits. 


Operational  Control 
Control  Surface 


Stabilator  (5c/ ^d) 
Aileron  (5^) 
Rudder  (5,) 


Deflection  Limits 
Positive  Limit 


15* 

20* 

30* 


(  8:26)  : 

Negative  Limit 


-29* 

-20* 

-30* 


Rotary  Balance  Database  limitations: 

Control  Surface  Positive  Limit  Negative  Limit 


Stabilator  (5c)  0*  -25* 
Differential 

Stabilator  (5d)  11*  -11* 
Aileron  (5^)  20*  -20* 
Rudder  (Sj.)  30*  0* 


Data  Analysis  and  Preparation 

In  examining  the  raw  data,  it  was  identified  that  the  rotary  balance  data  showed 
problems  in  the  static  configurations  (i.e.  Q.  =  0).  It  is  difficult  to  identify  what 
characteristics  of  the  data  is  due  to  poor  testing  or  is  inherent  in  the  character  of 
the  rotary  aerodynamics.  To  identify  possible  problems,  each  data  file  was 
compared  to  the  data  in  reference  (5)  to  identify  if  any  test  cases  showed  very 


15 


different  results.  All  coefficients  compared  fairly  well  except  for  a  few  isolated 
cases  of  and  Cy.  A  few  cases  showed  negative  signs  missing  on  a  series  of 
axial  force  coefficients  and  were  easily  corrected.  More  severe  problems  were 
encountered  with  particular  configuration  datasets  of  and  Cy.  The  magnitudes 
of  the  coefficient  data  for  configurations  1  and  4  on  Table  I  for  were  very 
different  than  cases  with  similar  configurations  in  both  reference  (5)  and  the 
available  data  on  floppy  disks.  These  cases  were  simply  eliminated  since  no  other 
data  was  available  to  replace  those  configurations.  Similarly  the  test  data 
configurations  1  and  13  for  Cy  also  showed  very  different  results.  However,  the 
data  for  case  1  in  reference  (5)  was  satisfactory  (compared  relative  to  the  1988 
Aerobase  coefficient  data)  and  was  used  as  a  replacement.  In  general,  the  overall 
estimations  of  the  axial  and  side  force  coefficients  are  questioned. 

As  the  experimental  data  was  collected  during  the  spin  tunnel  testing,  static 
results  were  obtained  prior  to  both  clockwise  and  counter  clockwise  rotations. 
Therefore,  twice  as  much  data  was  available  with  static  conditions  (Qss=0)  as  with 
any  other  rotation  rate.  This  provided  an  opportunity  to  examine  the  consistency 
of  the  coefficient  values  for  static  condition.  Many  test  cases  showed  reasonably 
close  coefficient  values  however  most  cases  were  inconsistent.  Having  no 
available  means  to  decipher  which  was  the  best  coefficient  value  from  the 
clockwise  and  counter  clockwise  tests,  the  static  condition  coefficient  data  were 
simply  averaged. 

In  preparation  for  creating  a  hybrid  model,  it  needed  to  be  determined  whether 
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the  two  coefficient  databases  were  reasonably  based  on  the  same  testing 
conditions  before  blending  them.  The  conditions  chosen  for  the  1988  F-15 
Aerobase  model  has  been  a  thrust  setting  of  8,300  lbs  (military  power,  thrust 
setting  for  trim  conditions,steady  level  flight)  at  0.6  Mach  and  20,000  feet.  The 
testing  conditions  for  the  rotary  balance  coefficients  were  obtained  in  a  low 
Reynolds  number  environment  which  was  not  specifically  set  for  variable  flight 
conditions.  The  CG  location  as  well  as  the  aircraft  body  fixed  coordinate  frames 
were  identical  for  both  databases.  Both  databases  are  body  frame  referenced  with 
the  (CG)  position  of  25.65  percent  mean  aerodynamic  chord.  Figure  2.1  and  2.2 
show  comparison  of  the  1988  F-1 5  Aerobase  as  modelled  by  Baumann  (7)  and  the 
raw  rotary  balance  force  and  moment  coefficient  data  for  a  selected  configuration. 
Of  the  six  coefficients  (three  force  and  three  moment),  the  moment  coefficients 
have  the  better  correlation  to  the  1988  Aerobase.  However,  given  the  1988 
Aerobase  and  rotary  balance  data  base  were  acquired  from  different  testing 
facilities,  have  possible  variation  in  Reynolds  number  and  possible  testing  errors 
may  be  evident  in  the  spin  tunnel  data,  the  coefficients  would  not  be  expected  to 
match.  Some  coefficients  as  with  the  axial  force  and  roll  moment  coefficients  show 
large  qualitative  differences.  However,  their  magnitude  differences  are  everywhere 
relatively  small.  Results  may  not  be  sensitive  to  the  differences.  The  effects  of 
variations  of  the  axial  and  side  force  coefficient  were  investigated. 

Data  analysis  is  a  highly  creative  venture.  Most  engineering  and  research  data 
are  assembled  in  unbalanced  sets  that  provide  no  attention  to  the  requirements 
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Figure  2.1  Comparison  of  force  coefficients  for  the  raw  rotary  balance  data 
and  the  1988  F-15  Aerobase  model  by  Baumann  (7). 
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Figure  2,2  Comparison  of  the  moment  coefficients  for  the  raw  rotary  balance 
data  and  the  1988  F-15  Aerobase  model  by  Baumann  (7). 
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of  statistical  analysis  of  the  experimental  results.  However,  because  of  the  limited 
availability  of  data,  the  analyst  must  work  with  what  he  gets.  The  bifurcation 
analysis  package  AUTO  (discussed  in  Chapter  III)  works  best  with  functions  with 
continuous  first  partial  derivatives.  Previous  works  using  AUTO  (6), (7), (8), (24) 
demonstrated  problems  performing  calculations  with  functions  that  have 
discontinuities.  It  was  determined  to  develop  a  polynomial  as  a  function  of  seven 
variables  for  each  coefficient  providing  a  representation  of  the  coefficients  without 
any  discontinuities  or  sharp  changes.  This  enabled  the  full  character  of  the 
database  to  be  modelled  in  terms  of  how  the  coefficients  varied  with  changes  of 
the  control  surfaces  as  well  as  angle-of-attack,  angle  of  sideslip  and  rotation  rate, 
Representing  the  data  in  terms  of  polynomial  equations  also  provided  an  efficient 
means  for  computation.  Table  look-ups  and  cubic  spline  interpolation  between 
discrete  values  were  considered.  However,  the  high  order  characteristics  of  the 
data  and  limited  depth  in  configurations  led  concerns  for  sudden  transitions  or 
discontinuities  in  the  coefficient  values  resulting  in  problems  with  AUTO.  As  well, 
there  are  extreme  complexities  of  implementing  a  seven  variable  local  fit  cubic 
spline  routine.  A  spline  routine  used  for  three  dimensional  imaging  was 
researched  and  examined  which  led  to  the  conclusion  that  a  seven  variable  effort 
would  be  too  complex  for  the  intention  of  this  research.  Polynomial  curve  surface 
fitting  of  the  data  will  introduce  smoothing  and  possibly  invite  characteristics  into 
the  data  that  are  not  discretely  present  in  the  raw  data.  However,  the  developed 
characteristics  will  be  driven  by  the  pattern  established  by  the  raw  data.  For  the 
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reader  interested  in  discussion  of  problems  and  methods  of  fitting  equations  to 
data,  reference  (12)  is  recommended. 

Curve  fitting  is  an  art.  With  the  size  of  the  database  (approximately  200Q 
experimental  cases),  the  odds  of  developing  the  polynomial  that  perfectly  fits  are 
very  small.  However  rewarding  results  could  be  obtained  even  with  a  moderate 
fit.  Any  software  used  should  provide  the  user  with  control  over  the  analysis; 
providing  the  user  with  a  variety  of  options.  The  need  for  a  curve  fitting  software 
that  could  handle  a  function  of  seven  variables  limited  the  choice  to  two  software 
packages,  SAS  (31 )  and  STATISX  (36).  SAS  is  a  mainframe  based  package  and 
STATISX  is  a  PC  based  package.  Both  software  were  used  for  the  development 
of  the  polynomials.  SAS  and  STATISX  were  able  to  provide  statistical  measures 
of  how  accurate  the  curve  fit  was.  They  also  provided  direct  comparison  between 
the  raw  data  and  the  polynomial  function  fit. 

The  raw  data  consisted  of  many  outliers  which  were  difficult  to  identify.  Plots 
of  various  dimensions  of  the  data  were  made  to  assist  in  outlier  identification  and 
to  identify  its  complexity  for  polynomial  fitting.  Even  though  most  datasets  had 
non-dimensional  rotation  rates  to  ±  0.7,  it  Was  determined  that  the  values  of  the 
coefficients  at  rates  above  ±  0.5  were  inconsistent  in  sign  and  magnitude  and 
would  contribute  more  difficulty  to  the  SAS  fitting  than  contribute  to  the  modelling 
of  the  database. 

The  limited  configurations  available  in  the  rotary  balance  database  determined 
the  maximum  order  of  any  of  the  variables  in  the  polynomial  model.  Table  III 
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indicates  the  maximum  order  of  each  of  the  variables  based  on  the  available 
variation  of  each.  It  was  detennined  to  limit  the  polynomial  dimension  to  fourth 
order  for  any  individual  or  coupled  variables.  Higher  order  polynomials  were 
possible  however  the  complexity  in  determining  a  quality  fit  was  increased.  The 
increase  in  data  processing  depth  seemed  to  outweigh  the  benefit  of  a  closer  fit. 
A  sixth  and  eighth  order  model  were  developed  for  some  of  the  coefficients.  The 
higher  order  models  did  not  provide  a  tremendous  improvement  to  the  fit  and 
displayed  unwanted  behavior  with  some  of  the  independent  variables.  The 
complexity  in  identifying  the  cause  of  the  behavior  outweighed  the  improvements. 

Table  III  Maximum  Order  of  Each  Polynomial  Variable 

Variable  Maximum  Polynomial  Order 

8 
2 

10 
2 
3 
1 
2 

The  curve  fitting  began  by  identifying  all  combinations  of  tne  variables  to  a 
fourth  order  multi-variable  polynomial.  Some  combinations  should  not  even  be 
considered.  Control  surface  effects  would  not  be  strongly  coupled  since  their 
positions  are  mechanically  independent  and  aerodynamic  coupled  influences  are 
negligible.  For  example,  aileron  settings  may  create  aerodynamic  influences  on 
rudder  derivatives  but  the  distance  between  the  control  surfaces  result  in  the 
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effects  being  negligible.  The  only  coupled  terms  involving  multiple  control 
deflections  are  the  differential  stabilator  and  aileron  since  they  are  inherently 
coupled  in  the  design  of  the  aircraft  by  the  linear  relation,  6d=(0.3)6a.  Six  of  the 
sixteen  available  configurations  in  Table  I  varied  independent  of  the  defined 
linear  relationship  with  83.  This  allowed  the  resulting  polynomial  to  have 
independent  terms  of  8^  in  addition  to  coupled  with  83.  It  was  left  to  SAS  to  identify 
terms  that  had  no  correlation  to  the  behavior  of  the  data.  Indicating  the  strength 
of  SAS  and  giving  validity  to  the  polynomial  "surface"  fit,  SAS  was  able  to  identify 
those  coupled  terms  that  were  not  needed  in  the  polynomial.  The  coupled  term 
of  aileron  and  differential  stabilator  deflection  was  identified,  as  predicted. 

STATISX  was  used  to  prepare  the  raw  data  for  use  in  SAS.  The  SAS  routine 
GLM  was  used  to  perform  a  linear  regression  to  determine  the  coefficients  of  the 
defined  polynomial.  Using  a  method  of  maximizing  the  standardized  residual  mean 
square  value  towards  1.0  and  identifying  strong  terms  using  the  student’s  t 
distribution  function,  the  "best"  polynomial  was  determined.  At  times,  the  process 
involved  trial  and  error  in  removing  and  adding  terms.  However,  each  SAS  run 
presented  evidence  of  the  sensitivity  of  a  coefficient  to  particular  variables  or  their 
combination.  The  student’s  t  distribution  function  provided  "clues"  to  what  the 
next  "guesses"  should  be.  The  residuals  of  each  polynomial  fit  were  plotted 
against  the  independent  variables  to  identify  trends  and  to  assist  in  determination 
of  the  "quality"  of  the  fits. 

As  an  experiment,  to  test  the  method  used,  a  fourth  order  polynomial  of  four 
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variables  was  created  and  a  random  database  of  results  was  constructed. 
Knowing  the  "answer"  to  the  polynomial  fit,  SAS  runs  were  performed  utilizing 
similar  decisions  made  on  the  rotary  balance  database.  This  was  performed  to 
validate  the  method  utilized.  The  resulting  polynomial  was  very  similar  to  the  one 
created  showing  validity  in  the  development  process  used  for  the  rotary  balance 
coefficient  polynomials. 

The  first  attempt  at  fitting  the  data  was  to  break  each  coefficient  into  6  individual 
polynomials  broken  at  20°  increments  with  10°  overlapping  in  angle  of  attack.  The 
objective  was  to  reduce  the  amount  of  data  each  fitting  needed  to  use  and 
decrease  the  residuals  for  a  more  accurate  representation  of  the  data.  The 
decision  to  use  angle  of  attack  was  determined  by  examining  the  effects  of 
breaking  up  the  data  base  into  sets  based  on  each  of  the  seven  independent 
variables  (i.e.  dividing  the  entire  database  into  three  groups  defined  by  angle  of 
sideslip:  p=-10°,p=+l0°,p=+10°).  Seven  combinations  of  the  data  (22  resulting 
datasets)  were  developed  representing  each  of  the  independent  variables.  Using 
the  complete  fourth  order  polynomial  representation  (approximately  200  terms),  the 
standardized  residual  mean  square  results  were  compared  to  see  which 
combination  of  datasets  provided  the  "best"  fit.  This  processing  also  assisted  in 
identifying  which  variables  were  most  difficult  to  model.  Datasets  divided  by  angle 
of  attack  showed  the  best  standardized  residual  mean  square  values  and  those 
by  angle  of  sideslip  the  worst.  This^resulted  in  the  development  of  36  equations 
representing  the  six  coefficients.  The  resulting  polynomials  consisted  of  equations 
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ranging  from  15  to  37  multi-variable  terms.  Once  the  36  equations  were 
developed,  79  distinct  monomials  in  7  variables  were  identified.  A  matrix  of  the 
polynomial  coefficients  was  constructed  with  those  terms  not  required  for  a 
particular  polynomial  term  given  a  null  value.  Each  set  of  six  equations  was 
blended  together  using  a  cubic  spline  relationship  similar  to  Eqn  (45)  of  Chapter 
IV. 

When  the  36  equation  model  was  used  with  AUTO,  each  of  the  six  coefficients 
encountered  problems  in  performing  parameter  sweeps.  This  indicated  difficulties 
with  the  transitions  between  the  individual  equations.  As  shown  in  Figure  2.1 , 
Cj^C^  and  C„  show  fairly  linear  behavior  which  resulted  .in  minimal  problems  with 
transitions  from  individual  equations  in  their  modelled  polynomials.  However,  each 
did  experience  occasional  problems.  Cy  (see  Figure  2.5), C,  and  had  much 
more  character  for  the  polynomials  to  represent  which  resulted  in  more  frequent 
discontinuous  transitions  between  the  equations  representing  the  coefficient. 
Figure  2.3  shows  the  C,  and  six  part  polynomial  models  for  two  different 
configurations  as  compared  to  the  raw  rotary  balance  data.  The  model  obviously 
was  an  excellent  fit.  However,  the  sharp  transitions  between  the  Individual 
equations  are  evident.  Figure  2.3  is  not  representative  of  the  "worst"  case 
discontinuities.  Because  each  equation  of  a  coefficient  polynomial  was  defined  on 
a  limited  database,  each  displayed  different  behavior  to  the  independent  variables. 
Most  often  the  behavior  difference  was  slight.  Occasionally  gross  differences 
occurred.  It  was  difficult  to  examine  every  dimension  of  the  seven  variable 
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polynomials  and  identify  when  difficult  transitions  occurred.  AUTO  encountered 
discontinuity  problems  on  every  continuation  run,  hindering  the  investigation.  Two 
and  one  half  months  of  development  were  invested  in  the  36  equation  model 
however  a  different  representation  of  the  data  had  to  be  pursued.  The  range  of 
angle  of  attack  for  blending  the  transition  between  equations  was  expanded 
however  the  results  with  AUTO  were  still  not  successful.  It  is  anticipated  that 
problems  with  individual  polynomials  for  particular  combinations  of  independent 
parameters  may  have  undesirable  behavior.  With  further  development,  the  six  part 
polynomials  could  work  and  would  be  the  better  polynomial  representation  of  the 
data  because  of  its  close  fit.  However  if  noise  were  present  in  the  data,  the 
polynomial  would  also  be  amplifying  its  contribution. 

A  two  piece  equation  was  also  developed  which  encountered  similar  problems. 
Figure  2.4  shows  the  C,  and  coefficients  for  two  different  configurations.  The 
second  case  for  each  coefficient  demonstrate  how  the  polynomial  can  show  very 
different  behavior  and  result  in  a  poor  fit  while  having  good  correlation  on  other 
configurations  as  shown  in  the  first  case.  Rather  than  experience  similar  problems 
as  with  the  six  part  model,  it  was  decided  to  focus  on  a  single  equation  fit; 
eliminating  the  discontinuity  problems,  allowing  the  full  database  to  be  used  and 
providing  better  control  over  the  overall  sensitivity  of  a  polynomial  to  individual 
parameters. 

A  single  equation  representation  of  each  coefficient  was  developed  utilizing  the 
79  identified  terms  of  the  36  equation  development.  The  standardized  residual 
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(a.cl):  p=0^Q33b/2V„=O,53=O^S,=O^8e=O^6,=O”  (b.c):  p=10“,63=-20".5,=-6\5,=30“ 


Figure  2.3  Examples  of  the  six  equation  polynomial  representation  of  the 
rotary  balance  data  compared  with  raw  data. 

(a,d):  p=0^^233b/2V,  =0,83=0“, 6,=0^53=0^5,=0“  (b,c):  p=1 0“,83=-20“,8,=-6“,8,=30“ 


Figure  2.4  Examples  of  the  two  equation  polynomial  representation  of  the 
rotary  balance  data  compared  with  raw  data. 
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square  criteria  ranged  from  0.80  to  0.98  with  being  the  worst  and  the  best 
coefficient  data  fit.  All  coefficients  except  had  a  value  greater  than  0.9.  The 
resulting  polynomials  were  continuous  resulting  in  smooth  operation  with  AUTO. 

The  coefficient  functions  were  plotted  for  different  sample  aircraft  configurations 
and  states  comparing  them  to  the  raw  rotary  balance  data  and  the  1 988  Aerobase 
to  validate  the  accuracy  of  the  fits.  Figures  2.5  and  2.6  are  an  indication  of 
comparison  of  the  single  79  term  polynomial,  raw  rotary  balance  data  and  the 
Baumann  model. 

It  is  evident  from  Figures  2.5  and  2.6  that  the  polynomial  curve  fitting  has 
introduced  characteristics  that  are  not  evident  in  the  raw  data  as  well  as  missing 
some  that  are  evident.  What  the  polynomial  fitting  has  done  is  to  model  the 
database  as  a  whole  rather  than  looking  strictly  at  a  snap  shot  picture  of  the  data 
that  would  be  found  with  a  table  look-up  routine.  The  polynomial  fit  has  basically 
smoothed  the  data  pacifying  outliers  and  smoothing  subtle  changes  in  the  data. 
Comparing  the  rotary  balance  coefficient  polynomials  to  the  1988  Aerobase 
polynomial  fitted  data  shows  some  similarity  in  their  character  but  obvious 
differences.  It  is  apparent  that  the  differential  between  the  rotary  balance  data 
polynomial  fit  and  the  raw  rotary  balance  data  is  much  smaller  than  that  between 
the  rotary  and  1988  Aerobase  data.  During  spin  tunnel  testing,  the  problem  of 
aerodynamic  interference  caused  by  the  aircraft  model’s  wings/fuselage  vortices 
being  influenced  by  the  support  structure  can  be  significant.  In  addition,  any 
motion  that  alters  the  rotation  center  location  of  a  test  could  cause  error  in  the  post 
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Figure  2.5  Comparison  of  the  force  coefficient  rotary  balance  data  polynomial 
fitted  data,  raw  RB  data  and  the  Baumann  model  (7). 
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Figure  2.6  Comparison  of  the  moment  coefficient  rotary  balance  polynomial 
fitted  data,  raw  RB  data  and  the  Baumann  Model  (7). 
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processing  which  was  based  on  a  predetermined  rotation  center  location  (2:7-1 1). 
The  behavior  of  the  rotary  balance  data  in  Figures  2.5  (b)  and  2.6  (a)  could  be 
attributed  to  these  possible  inaccuracies  in  spin  tunnel  testing.  It  is  anticipated  that 
the  results  using  either  database  exclusively  or  combined  will  provide  quite 
different  results.  The  Baumann  (7)  1988  Aerobase  coefficient  representation  and 
the  79  term  rotary  balance  database  polynomial  can  be  referenced  on  page  119 
and  145  respectively  in  Appendix  C. 
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III.  Bifurcation  Theory 


This  chapter  will  discuss  the  basic  principles  of  bifurcation  theory  and  its 
application  toward  nonlinear  systems.  The  qualitative  study  of  nonlinear  differential 
equations  is  concerned  with  howto  deduce  important  characteristics  of  the  solution 
without  actually  solving  the  equations.  Through  bifurcation  theory  this  study  can 
investigate  highly  nonlinear  motion  of  aircraft  in  a  spin  without  the  need  to  basically 
destroy  information  contained  in  the  nonlinear  equations  by  linearizing  them  with 
small  perturbation  analysis.  There  are  many  types  of  bifurcations  and  each  type 
has  a  different  effect  on  the  response  of  a  system.  The  concepts  to  be  discussed 
are:  equilibrium  points,  phase  space,  stability,  turning  points,  bifurcation  points, 
periodic  solutions,  Hopf  points  and  an  example  of  unfolding  of  an  organized  center. 
Most  of  the  information  in  this  chapter  is  referenced  from  Seydel  (34).  Other 
useful  texts  on  the  subject  of  bifurcation  theory  are  references  (15)  and  (21).  A 
brief  description  of  the  software  program  AUTO  will  follow  with  a  short  discussion 
on  homotopy,  an  important  application  of  continuation  theory  needed  for  this 
research. 

Equilibrium  Points 

When  a  system  is  in  physical  equilibrium  or  in  a  steady  state  for  bodies  in 
motion  (  i.e.  an  aircraft),  the  states  that  describe  the  system  are  termed  the 
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equilibrium  points.  Equilibrium  points  are  also  referred  to  as  stationary  points.  In 
the  study  of  aircraft  dynamics,  only  autonomous  systems  will  be  analyzed.  For  an 
autonomous  system,  the  differential  equations  of  motion  do  not  explicitly  contain 
time  on  the  right  hand  side  of  the  equations.  An  autonomous  system  can  be 
written  as 

u  =  ^u)  0) 

where  u  is  an  n-dimensional  vector  of  state  variables  and  f  is  an  n-dimensional 
vector  of  functions  describing  the  motion  of  the  system.  The  system  is  in 
equilibrium  when  the  states  are  constant,  u  =  0.  The  states  that  describe  this 
equilibrium  of  the  system  would  satisfy  the  equation 

f{u)  =  0  (2) 

These  states  are  called  the  equilibrium  points  or  stationary  points. 

Phase  Space 

Suppose  the  state  of  a  system  is  described  by  the  state  vector  u  and  the 
nonlinear  equation 

u  =  fiu,},)  (3) 

defining  the  behavior  of  the  system  where  a.  is  a  control  parameter  of  the  system. 
Points  along  the  solution  of  Eqn  (3)  consist  of  a  time  coordinate  t,  the  fixed  value 
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control  parameter  X  and  the  n  dimensional  space  coordinates  (yi,y2,—.yn)-  the 
vector  function  f  is  continuously  differentiable,  there  is  a  unique  trajectory  through 
each  point  Uq  of  Eqn  (3). 

Consider  a  three  dimensional  phase  space  representing  the  coordinates 

(yi.y2.^)-  At  equilibrium,  X=0  ,  for  a  fixed  value  of  X  the  solution  will  remain  in  a 

plane,  the  phase  plane.  A  stationary  solution  is  represented  by  a  point.  Figure  3.1 . 
As  the  parameter  X  is  varied,  the  position  of  the  state  solutions  change  in  the 
phase  plane.  If  the  parameter  A.  is  freely  variable,  the  equilibrium  solution  path 
forms  a  curve  in  the  (yi.yj.A.)  space.  Figure  3.2.  These  curves  represent 
equilibrium  solutions  as  they  depend  on  variable  parameters  of  the  system.  This 
is  simple  to  visualize  in  three  dimensions  but  consider  a  set  of  eighth  state 
equations  of  motion  of  an  aircraft  dynamic  model  and  the  multiple  parameters  that 
can  be  varied.  The  parameter  X  could  represent  the  control  deflections  of  the 
aircraft’s  rudder,  stabilator,  ailerons  or  any  parameter  within  an  aircraft  attitude 
control  system.  It  becomes  evident  that  the  visualization  of  the  equilibrium  solution 
curves  can  become  quite  complex.  However,  qualitative  information  can  be  found 
by  viewing  the  curve  projections  of  a  chosen  system  state  as  the  value  of  a  single 
control  parameter  is  varied,  i.e.  viewing  the  aircraft  state  a  as  elevator  deflection 
is  varied.  The  resulting  projection  is  a  bifurcation  diagram.  Applications  of  a 
bifurcation  diagram  utilizedTor  this  investigation  are  the  examination  of  effects  on 
particular  states  of  aircraft  dynamics  as  control  surfaces  are  varied.  The  resulting 
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Figure  3.1  A  stationary  solution  in  several  phase  planes  (41 :3). 


Figure  3.2  Solution  curves  in  (y,X)-  space  (41:3). 
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solution  curves  provide  information  for  identification  of  effective  controls  for  aircraft 
maneuvers  or  of  unusual  aircraft  states  that  could  occur  when  the  aircraft  is 
configured  with  particular  control  surface  deflections.  Figure  3.3  is  an  example  of 
a  bifurcation  diagram.  Bifurcation  diagrams  also  allow  for  the  identification  of 
multiple  states  a  system  can  attain  for  a  given  value  of  a  control  parameter,  X.  It 
is  advantageous  to  identify  configurations  that  could  lead  to  abrupt  changes  in  the 
aircraft  dynamics  as  the  aircraft  states  jump  to  a  different  solution  value. 
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Figure  3.3  Example  Bifurcation  Diagram. 


Stability 

Stability  is  the  tendency  of  a  system,  when  disturbed  from  a  given  equilibrium, 
to  return  to  that  equilibrium.  An  equilibrium  may  be  stable  for  a  small  perturbation 
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but  unstable  for  a  large  perturbation  showing  that  a  system’s  stability  is  only  valid 
within  the  domain  of  attraction  of  the  solution  of  interest.  Figure  3.4  is  a  way  of 
looking  at  the  various  definitions  of  stability.  A  marble  placed  perfectly  balanced 
on  the  top  of  the  hill  (A)  is  an  unstable  equilibrium  state.  Any  disturbance  will 
cause  the  marble  to  leave  its  position  and  never  return.  At  position  (B)  the  marble 
is  apparently  stable  for  small  disturbances  however  under  a  strong  influence,  the 
marble  may  leave  this  location.  Equilibrium  locations  (B)  and  (C)  represent  weak 
stable  equilibrium  and  location  (D)  represents  a  strong  stable  equilibrium.  A  larger 
disturbance  at  location  (D)  will  allow  return  of  the  marble  to  its  original  stable 
position. 


Stability  of  a  stationary  point  determines  whether  the  state  of  the  system  is 
attracted  to  or  repelled  from  the  point.  Local  stability  refers  to  the  stability  in  a 
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small  region  around  a  stationary  point.  Local  stability  of  a  fixed  point  of  a 
nonlinear  system  can  be  calculated  by  determining  the  eigenvalues  of  the 
equivalent  linear  system.  The  eigenvalues  are  obtained  by  linearizing  the 
equations  of  the  system  about  the  stationary  point  of  interest.  A  stationary  point 
is  stable  if  the  real  parts  of  all  eigenvalues  are  less  than  zero  and  unstable  if  any 
eigenvalue  has  a  real  part  that  is  greater  than  zero.  In  two  dimensions,  Figure  3.5 
shows  the  various  stability  types  in  the  phase  space.  If  the  eigenvalues  are  real 
numbers  greater  than  zero,  the  stationary  point  is  a  source  or  if  less  than  zero  a 
sink.  If  the  eigenvalues  are  real  numbers  with  opposite  sign,  a  saddle  point  is 
formed.  For  complex  conjugate  eigenvalues  with  nonzero  real  parts,  if  Real  <  0 
a  stable  spiral  is  formed  else  if  Real  >  0  an  unstable  spiral.  For  the  condition  with 
complex  conjugate  eigenvalues  with  zero  real  parts  a  center  is  formed.  With  multi- 
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dimensional  systems,  the  higher  dimensional,  analogous  phase  behavior  to  Fig  3.5 
must  be  imagined.  The  behavior  are  too  complex  to  visualize  however  do  occur. 

Turning  Point 

A  simple  example  of  a  turning  point  can  be  introduced  using  the  scalar  equation 

y  =  A.  -  (^) 

The  equilibrium  solutions  can  be  determined  from 

0  =  A.  -  y2  (5) 

The  solutions  y(A.)  of  equation  (4)  form  a  parabola  that  is  only  defined  for  A>0.  At 
A.=0  there  is  only  one  solution  (y=0)  where  as  for  A>0  there  are  two  solutions.  The 
point  where  a  solution  begins  to  exist  {A.=0,y=0)  is  a  turning  point  or  also  referred 
to  as  a  limit  point.  Figure  3.6  graphically  shows  a  limit  point.  A  solid  line  depicts 
a  stable  solution  branch  and  a  dashed  line  an  unstable  branch.  This  convention 
will  be  used  throughout  this  report.  It  should  be  noted  that  in  n-dimensional 
systems  a  turning  point  does  not  always  separate  stable  equilibrium  from  unstable 
equilibrium.  An  eigenvalue  always  changes  sign  at  the  turning  point  but  others 
may  already  be  greater  than  zero. 

Turning  points  often  arise  in  pairs  resulting  in  hysteresis  effects.  Figure  3.7  is 
an  example  of  a  hysteresis  effect.  Characteristic  of  hysteresis  are  jump 
phenomena  which  take  place  at  A.,  and  A^.  As  A.  is  increased  along  the  upper 
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Figure  3.6  Representation  of  a  limit  point. 


Figure  3.7  Example  of  a  hysteresis  or  jump  phenomenon. 
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branch  approaching  the  right  limit  point,  increasing  the  value  of  X  beyond  the  value 
X2  will  cause  a  jump  to  the  lower  branch.  Similar  characteristics  can  be  found  for 
the  value  X,.  This  phenomenon  is  also  referred  to  as  a  jump  phenomenon  Such 
behavior  has  been  found  in  aircraft  dynamics  as  discussed  in  reference  (33)  by 
Schy  and  Hannah. 

Bifurcation  Points 

A  bifurcation  is  a  point  where  something  is  divided  into  two  parts.  A  bifurcation 
occurs  in  a  system  when  the  variation  of  an  independent.control  parameter  creates 
a  point  where  the  behavior  of  the  system  can  assume  one  of  two  different  states 
for  the  same  set  of  system  parameters.  Consider  a  system  described  by  the 
scalar  equation 

y  =  Xy  -  (®) 

Foi  all  values  of  X  there  is  the  trivial  solution  y=0.  If  ^5*0  the  non-trivial  solution  is 
y=X.  A  stability  exchange  occurs  as  depicted  in  Figure  3.8.  The  point  at  X=0  is 
referred  to  as  a  transcritical  bifurcation  point.  Other  qualitative  types  of  simple 
bifurcations  are  possible  but  will  not  be  discussed. 

When  the  branch  y=0  loses  stability  at  the  bifurcation  point  and  branches  into 
two  stable  trajectories,  the  event  is  referred  to  as  a  supercritical  pitchfork.  If  the 


39 


Figure  3.8  Example  transcritical  bifurcation  point. 

branch  y=0  gains  stability  at  the  bifurcation  point  and  the  bifurcation  of  unstable 
branches  occur,  the  result  is  called  a  subcritical  pitchfork.  Examples  of 
supercritical  and  subcritical  pitchforks  are  shown  in  Figure  3.9. 


Figure  3.9  Examples  of  supercritical,  subcritical  bifurcations. 
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Limit  Cycle 


A  limit  cycle  represents  regular  motion  such  as  the  vibration  of  a  string  or 
current  flowing  through  an  electrical  circuit.  A  limit  cycle  is  an  isolated  periodic 
solution  of  an  autonomous  system,  represented  in  the  phase  plane  or  phase  space 
for  n-dimensional  systems  by  an  isolated  closed  path.  Unlike  a  center,  the 
neighboring  paths  are  not  closed  but  either  spiral  into  or  away  from  the  limit  cycle. 
Figure  3.10  is  an  illustration  of  a  stable  limit  cycle.  For  the  stable  limit  cycle,  any 
state  of  the  system  near  the  limit  cycle  will  drift  into  the  periodic  motion  defined  by 
the  limit  cycle. 


Figure  3.10  Example  stable  limit  cycle. 


Hopf  Bifurcation 

Hopf  bifurcation  is  the  door  that  opens  from  the  small  room  of 
equilibria  to  the  large  hall  of  periodic  solutions  (34:61). 
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The  type  of  bifurcation  that  connects  equilibria  with  periodic  motion  is  Hopf 
bifurcation.  Consider  the  system 


y,  =  ^yt  *  Vi  -  yi(yi^  yl) 

%  =  -y,  +  'ky^  -  yaCy,  yg) 

y^=y2=0  are  the  only  equilibrium  values  for  all  k  and  no  stationary  bifurcations 
occur.  The  Jacobian  matrix  for  this  system  is 


k  -1 
1  k 


and  has  eigenvalues  k=£\.  The  equilibrium  solutions  are  stable  for  A.<0  and 
unstable  for  X>0  with  a  loss  of  stability  at  X=0.  Using  polar  coordinates 


y,  =  rcosG  ,  =  /sin0 

the  system  of  equations  can  be  simply  expressed  as 

fi) 

and  0  =  1 


(9) 


(10) 


If  X<0  the  entire  phase  diagram  is  a  stable  spiral.  If  X>0  an  unstable  spiral  is 
formed  at  the  origin  surrounded  by  a  stable  limit  cycle  which  grows  out  from  the 
origin  as  k  increases.  The  origin  changes  from  being  asymptotically  stable  to 
being  unstable  without  passing  through  the  stage  of  a  center.  Figure  3.1 1  shows 
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the  phenomenon  in  two  dimensional  slices  and  Figure  3.12  in  three  dimensions  by 
including  X.  On  the  bifurcation  diagram,  periodic  solutions  will  be  depicted  by 
circles  along.a  branch  as  shown  in  Figure  5.2(c)  of  Chapter  V.  Closed  circles  will 
depict  stable  limit  cycles  and  open  circles  unstable  limit  cycles.  It  should  be 
understood  that  a  stable  periodic  orbit  is  approached  by  nearby  trajectories, 
whereas  trajectories  leave  a  neighborhood  of  an  unstable  periodic  orbit. 


Figure  3.11  Development  of  a  limit  cycle  in  a  Hopf  bifurcation. 


Figure  3.12  Three-dimensional  representation  of  a  Hopf  bifurcation. 
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Organized  Centers 


Locations  in  solution  branches  where  dynamical  behavior  changes  in  the 
sequence;  no  multiplicity,  then  two  jumps,  then  four  jumps,  then  two  jumps,  is  one 
example  of  an  organized  center.  This  phenomena  can  best  be  explained 
graphically.  Figure  3.13  shows  various  bifurcation  diagrams  as  a  second  control 
parameter  7  is  varied.  For  this  example  Y5>74>Y3>Y2>Yi-  Fo''  Ta  ^wo  branches  are 
found  without  connection,  the  upper  branch  is  an  isolated  branch,  the  path  A-B  is 
smooth.  For  73  the  situation  is  the  same  but  the  two  branches  are  closer.  For  74 
there  is  no  longer  an  isolated  branch,  with  the  overall  branch  structure  now 
resembling  a  mushroom.  In  this  state,  the  consequences  for  the  path  from  A  to 
B  is  severe.  There  are  two  hysteresis  jumps.  There  is  a  transcritical  bifurcation 
that  occurs  between  the  values  73  and  74  that  separates  a  two  jump  situation  to  a 
four  jump  situation.  At  some  value  of  7  the  branches  merge  in  an  isola  point.  An 
isola  point  would  be  equivalent  to  the  transcritical  bifurcation  point  in  Fig  3.8  with 
two  of  the  like  stability  branches  being  a  closed  path.  Increasing  7  further  can 
result  in  a  two  jump  situation.  This  representative  sequence  of  high-order 
bifurcations  is  governed  by  theory  beyond  the  scope  of  this  investigation.  For 
discussion  on  understanding  why  these  events  occur  and  further  examples  of  this 
phenomena,  Seydel  (34)  is  referred.  This  type  of  phenomena  occurs  in  reality  and 
has  been  shown  during  this  investigation. 
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AUTO  Software 


The  software  tool  utilized  in  this  investigation  for  continuation  and  bifurcation 
analysis  is  AUTO.  AUTO,  written  by  Eusebius  Doedel  of  Concordia  University,  is 
a  collection  of  FORTRAN  routines  concerning  the  nurherical  analysis  of  nonlinear 
systems  of  algebraic  and  ordinary  differentiaf  equations.  The  software’s  primary 
purpose  is  to  compute  branches  of  stable  or  unstable  periodic  solutions  of  systems 
of  the  form  of  Eqn  (3).  Given  a  function,  the  Jacobian  of  the  function,  the 
derivative  f;^,  a  steady  state  solution  for  some  value  X  and  a  number  of  control 
parameters,  AUTO  can  compute  steady  state  branches,  accurately  determine 
steady  state  and  Hopf  bifurcation  points  and  switch  branches  at  such  points.  The 
key  to  initiating  an  analysis  using  AUTO  is  the  identification  of  a  steady  state 
solution  of  the  system  of  equations  to  be  analyzed.  This  can  be  accomplished  by 
considering  a  simple  flight  condition  or  the  methods  of  homotopy  could  be 
exercised.  Described  is  only  a  partial  outline  of  AUTO’s  capabilities.  The  reader 
is  encouraged  to  reference  the  AUTO  User’s  Manual  (13)  for  additional  information 
on  its  capabilities  and  application. 

Continuation  and  Homotopy 

Continuation  theory  is  the  methodology  of  answering  the  question  of  how 
solutions  of  equations  vary  with  a  parameter.  Continuation  methods  involve  four 
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basic  elements:  (1)  an  initial  guess  of  the  solution,  (2)  choice  of  the  system 
parameter  to  vary,  (3)  a  correction  iteration  of  the  solution  and  (4)  variation  of  the 
iteration  step  size^  Given  a  set  of  differential  equations  describing  a  dynamical 
system,  a  numerical  integration  routine  such  as  Euler’s  or  Newton’s  method  could 
be  used  to  integrate  the  equations.  Predictions  of  the  next  solutions  along  a 
branch  are  determined  as  small  incremental  changes  are  taken  of  the  chosen 
system  parameter.  A  corrector  iteration  is  then  used  improving  the  guess  of  the 
next  solution.  As  the  solution  branches  are  identified,  variable  step  lengths  are 
used  to  allow  for  small  details  to  be  resolved  and  not  skipped,  assist  in 
transitioning  through  turning  points  and  to  aid  in  processing  of  the  numerical 
iteration  methods.  The  result  is  a  tracing  of  equilibrium  solutions  of  the  system 
of  differential  equations. 

One  important  application  of  continuation  is  homotopy.  Consider  an  equation 
f(u)=0  which  is  difficult  to  solve  for  a  solution.  An  initial  guess  may  be  even  harder 
to  determine.  Note  that  this  equation  is  probably  nonlinear  and  would  require 
iterating  to  solve  .  Iterative  solution  algorithms  usually  converge  very  slowly  or 
diverge  away  from  the  solution  if  the  initial  guess  is  way  off.  Assume  that  an 
equation  g(u)=0  is  known  that  is  easily  solved  with  solution  Uq  and  can  be  obtained 
by  simplifying  f(u)=0.  Homotopy  is  a  construction  of  equations  that  are  linked 
together  and  are  solved  one  at  a  time.  The  last  equation  solved  is  the  original 
equation  f(u)=0.  The  solution  of  each  successive  equation  is  used  as  an  initial 
guess  for  the  next  equation.  This  describes  a  discrete  homotopy  with  finite 
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number  of  equations.  A  continuous  homotopy  is  where  a  parameter  A.  is 
introduced  that  varies  on  the  interval  0  <A.<1.  This  leads  to 

f[u,X)  =  0  ,  0  <  X  <  1 
f{u,0)  =  giu)  =  f{u) 

An  example  of  such  a  process  is  with 

fiu,l)  =  Xf{u)  +  (1  -  X)g{u) 

AtX=0,  Ug  is  the  solution.  As  X  is  varied  to  1.0,  the  solution  set  transitions  to  that 
of  the  equation  of  interest  f(u)=0.  The  method  of  homotopy  will  be  applied  in 
chapter  IV  to  acquire  equilibrium  solutions  required  by  AUTO  for  the  rotary  balance 
data  model  using  the  known  solutions  of  the  McDonnell  model  (24).  For  the 
reader  interested  in  a  short  tutorial  on  the  methods  of  continuation,  the  article 
"Tutorial  on  Continuation"  by  Seydel  (35)  is  highly  recommended. 


(11) 


(12) 
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IV.  Aircraft  Model 

The  aircraft  for  this  investigation  's  the  McDonnell  Douglas  Corporation  F-1 5B, 
a  two  place  high  performance  fighter  developed  for  the  United  States  Air  Force. 
The  F-1 5  was  designed  with  a  cockpit  sized  for  two  crew  members  requiring  the 
only  major  change  of  a  larger  canopy  to  form  the  two  seat  (F-15B)  from  a  single 
seat  version  (F-15A).  This  design  feature  allowed  for  application  of  the  F-15  in  a 
multitude  of  roles  in  training,  air-to-air  combat,  air-to-ground  combat  and 
reconnaissance.  The  flight  characteristics  of  the  F-1 5A  and  F-1 5B  are  very  similar 
even  with  the  added  weight  of  the  second  crew  member  configuration.  This  will 
enable  some  flight  test  data  of  the  F-1 5A  to  be  considered  in  the  analysis  of  the 
F-1 5B.  One  design  feature  sought  for  in  the  F-1 5  was  an  aircraft  capable  of  high- 
energy  maneuverability  in  turning,  accelerating  and  climbing  in  order  to  gain  a 
tactical  advantage  in  combat.  The  low  wing  loading,  very  high  thrust  design  of  the 
aircraft  adds  to  its  air  superiority  capability.  This  air  superiority  was  demonstrated 
in  the  F-15’s  multiple  roles  in  Operation  Desert  Storm.  A  figure  of  the  aircraft  with 
control  surface  sign  conventions  is  located  in  Appendix  B  along  with  aircraft 
physical  dimensions  and  specifications  in  Appendix  A. 

The  F-15B  has  multiple  independent  control  surfaces,  left  and  right  aileron,  left 
and  right  rudder,  left  and  right  stabilator,  speed  brake  and  variable  inlet  ramps. 
The  ailerons  and  stabilator  were  set  by  the  manufacturer  to  act  differentially  with 
the  linear  relationship  6d=(0.3)8a .  The  speedbrake  was  not  modeled  because  it 
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is  not  a  nominal  control  surface  used  during  high  angle  of  attack  maneuvers  and 
is  designed  for  retraction  at  angles  of  attack  above  15°.  The  aircraft  is  an 
inherently  stable  design  without  the  use  of  the  Control  Augmentation  System 
(CAS).  This  stability  characteristic  enabled  this  analysis  to  be  performed  without 
CAS  engaged.  The  aircraft  was  modeled  for  gear  up,  retracted  flaps  and  stores 
carriage.  The  Baumann/McDonnell  Model  aerodynamic  coefficients  were  modeled 
for  a  flight  condition  of  20,000  feet  and  a  Mach  number  M<0.6  (7).  The  testing 
conditions  for  the  rotary  aerodynamic  coefficients  fall  within  this  criterion. 

Model  Development 

The  development  of  the  aerodynamic  force  and  moment  coefficients  for  the 
F-15  model  used  by  McDonnell  (24)  are  documented  in  the  thesis  by  Baumann 
(7:20-21).  Baumann  curve  fitted  the  1988  F-15  Aerobase  which  represents  the 
F-1 5  aerodynamic  coefficients  for  static  and  forced  oscillation  testing  as  well  as  the 
inclusion  of  some  actual  flight  test  data.  Results  of  those  curve  fits  are 
documented  in  the  AUTO  Driver  program  presented  in  Appendix  C.  The  thrust 
contributions  to  the  force  and  moment  coefficients  used  in  this  research  are  those 
developed  by  McDonnell  (24:24).  The  feature  of  thrust  vectoring  will  be  included 
in  the  rotary  aerodynamic  model  however  its  effectiveness  will  not  be  pursued  in 
this  investigation.  The  rotary  aerodynamic  coefficients  obtained  by  the  NASA 
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Langley  Spin  Tunnel  are  represented  as 


=  Axial  Force!  {q  S) 

(13) 

Cy  =  Side  Force/iq  S) 

(14) 

=  Normal  Force!  {q  S) 

(15) 

C,  =  Rolling  Moment! (q  S  b) 

(16) 

=  Pitching  Moment! {q  S  c) 

(17) 

C„  =  Yawing  moment!  {q  S  b) 

(18) 

These  coefficients  are  used  to  represent  the  rotary  contributions  for  the 

development  of  a  hybrid  model.  The  coefficients  are  calculated  using  the 

measured  forces  and  moments  generated  on  the  test  model  during  rotary  balance 

testing.  To  include  thrust  contributions  for  the  stand  alone  rotary  balance  data 

model,  the  thrust  terms  developed  by  McDonnell  (24)  are  combined  with  the  rotary 

balance  coefficient  terms.  The  resulting  equations  are: 

=  {Axial  force  +  TJ!(q  S) 

(19) 

Cy  =  {Side  Force  +  Ty)!{q  S) 

(20) 

=  {Normal  force  +  T^)f{q  S) 

(21) 
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C,  =  {Rolling  moment  +  T^^djy  -  T^,dry)/{q  S  b) 

=  (Pitching  moment  +  T^d^^  -  '^zdf((ll(q  S  c) 

C„  =  {Yawing  moment  +  Tydj^  -  T^,djy  -  T^fijy)l{q  S  b) 


(22) 


(23) 

(24) 


The  development  of  the  equations  of  motion  for  an  airplane  are  outlined  by 
McRuer,  et  al  (25)  which  represent  a  conventional  set  of  differential  equations 
defining  the  dynamics  of  an  aircraft.  With  the  assumptions  of  a  rigid  airframe, 
inertial  earth  fixed  reference,  constant  mass  and  mass  distribution  of  the  aircraft 
and  constant  gravity,  the  resulting  equations  of  motion  are  a  ninth  order  set  of 
differential  equations  for  a  body  fixed  frarhe  of  reference.  The  aircraft’s  state  of 
motion  can  be  described  by  the  nine  state  variables  (a,  p,  p,  q,  r,  0,  (|),  y\r,  V„)  and 
the  deflections  of  any  control  surfaces  defining  the  moments  and  forces  acting  on 
the  aircraft.  Following  the  development  by  McRuer,  et  al  (30)  the  following 
equations  are  formed. 

Translational  acceleration  equations: 


B  =  -  cos9  sinB  +  r  cosa 

mV,  V„ 


+  cos6  sind)  cosB 

mV.,  ^  V„ 


C+-£-  cos9  cosd)  sinB  -  p  sina 
mV  V 

III 


V,  =  V.  -5 —  C  -  sine  cosa  cosp 

+  _5_£  C  +  cose  sind)  sinB 

^  K  J 


*  lAc  ^9 


cose  cos(j)  sina  cosp 


Rotational  acceleration  equations: 


+  RH  c,-^kc„ 
4 


1  -  i _ 'l  k  pq 

I  I 

Z  X 


l"xz 
1  -  -Ji 
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q  S  b 


C,^Cn 


/2„ 

• 

1  - 
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J  „ 

X  2^ 

Aircraft  orientation  relative  to  the  earth  inertial  reference  frame  (Euler  angles): 

0  =  g  cos(|)  -  rsincj) 

^  =  p  +  (  g,  sincj)  +  r  cos  (|) )  tan0 

V  =  (  g  sin(|)  +  r  cos(j) )  sec0 

By  the  definition  of  the  Euler  angles,  the  yaw  angle  y,  is  decoupled  from  the  rest 
of  the  equations  of  motion.  For  application  to  aircraft  dynamics,  the  usual 
convention  for  Euler  angle  rotations  is  the  sequence:  yaw,  pitch,  roll.  Because  the 
yaw  rotation  does  not  change  the  gravity  vector  relation  to  the  body  frame  z  axis, 
the  aircraft  can  be  modeled  without  the  ij;  orientation  equation  resulting  in  an 
eighth  order  model. 

A  more  accurate  model  of  the  dynamics  of  an  aircraft  would  be  to  include 
variation  in  thrust  due  to  the  motion  of  the  aircraft.  As  an  aircraft  enters  high- 


(31) 

(32) 

(33) 
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energy  dynamic  attitudes  the  state  of  thP  flow  conditions  about  the  air  intakes  of 
the  engines  change  resulting  in  the  effective  thrust  varying.  The  current  model  has 
the  aircraft  thrust  set  at  a  fixed  setting  of  8300  Ibf  for  an  altitude  of  20,000  ft.  The 
effects  of  variable  thrust  are  evident  in  Figures  (5-8)  and  (5-9)  of  reference  (29) 
showing  changes  in  the  equilibrium  solution  branches  and  the  effects  are  further 
discussed  in  reference  (30).  This  inaccuracy  in  the  current  model  is  noted 
however  for  comparison  to  previous  studies,  this  modification  will  not  be  pursued 
for  this  investigation. 

Total  Rotation  Vector 


When  incorporating  rotary  balance  data  into  a  model  it  should  be  referenced 
under  the  same  conditions  upon  which  it  was  collected.  During  rotary  balance  spin 
tunnel  testing,  the  model  is  rotated  at  a  constant  rate  about  an  axis  parallel  to  the 
free  stream  velocity  vector  of  the  tunnel.  The  test  model  sees  a  constant 
configuration  of  a,  P  and  free  stream  velocity.  Therefore,  the  rotary  balance  data 
should  be  referenced  using  the  component  of  the  total  rotational  rate  that  is  along 
the  free  stream  velocity  vector.  The  total  rotational  rate  is  defined  in  terms  of  the 
body  axes  rates  as 

Q.  =  p  b,  +  q  L  +  r  L 

(34) 

and  I  I  =  \/p^  + 
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For  spin  analysis,  the  total  rotational  vector  Q  is  projected  into  four  rotation  rates, 

the  steady  state  component  of  the  rotation  vector  along  the  free  stream  velocity 
vector  and  posc.Qosc>fosc  which  represent  the  body  frame  components  of  the  residual 
rotational  vector.  Referring  to  Figure  4.1,  these  components  are  defined  as 

n,.,  =  (  5  •  I/, )  V,  =  n„  v,  (35) 

Q.  =  Q.  -  Q.  .=  D  b.  +  Q  br.  +  r  b,  (3®) 

osc  rot  r^osc  ^osc  ^2  osc  ^2 


is  the  projection  of  the  total  rotation  vector  onto  the  velocity  vector  and 

represents  the  rotary  aspects  of  the  total  rotation  vector.  is  the  orthogonal 
component  of  the  total  rotation  vector  and  is  representative  of  the  oscillatory 
dynamics  inherent  in  the  total  rotation  vector.  By  the  process  of  vector  projections 
the  oscillatory  components  are  defined  as 


=  p  -  QjjCOsacosp 

(37) 

=  p  - 

(38) 

=  r  -  QjjSinacosp 

(39) 

Since  the  rotary  balance  data  was  experimentally  obtained  as  function  of  a,  P 
and  assb/2V„  and  additionally  through  development  of  the  polynomials  as  a 
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■  —  -  .  I 

Figure  4.1  Representation  of  the  rotation  vector  in  the  aircraft  body  frame  and 
wind  frame. 


Figure  4.2  Orientation  of  the  rotation  vector  when  oscillatory  contributions  are 
negligible. 
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function  of  8a,5d,5e  and  6„  the  key  to  the  proper  use  of  rotary  aerodynamic  data 
rely  on  the  appropriate  definition  of  the  steady  state  rotation*  rate  Hawkins 
(17:184)  used  Qjs  =  (pcosa  +  rsina)  which  restricted  his  analysis  to  zero  sideslip 
angles  or  the  assumption  of  very  small  perturbations  in  sideslip.  A  more  robust 
representation  is 

Q.^  =  (pcosa  +  /sina)cosP  +  qsinp  (40) 

V 

which  is  found  from  . 

It  is  evident  that  when  a  stand  alone  rotary  aerodynamic  model  is  used  it  is 

assumed  that  the  oscillatory  components  of  Q.  are  relatively  small  (posc.Qosc>fosc  ~ 

0)  compared  to  (=^ifot)  indicating  that  the  rotation  vector  is  closely  aligned  with 
the  free  stream  velocity  vector  (Figure  (4.2)).  If  the  orientation  of  the  total  rotation 
vector  were  grossly  misaligned  with  the  free  stream  velocity  vector,  the  oscillatory 
terms  would  play  a  larger  role  in  the  representation  of  the  dynamics  of  the  aircraft. 
This  would  support  the  argument  for  the  need  of  both  forced  oscillation  and  rotary 
aerodynamic  data  to  properly  represent  the  aircraft’s  motion.  However,  for  spins 
which  are  the  basis  for  this  investigation,  it  is  reasonable  to  assume  the  steady 
state  rotation  and  free  stream  velocity  vector  are  closely  aligned. 
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Rotary  Balance  Model 


For  the  stand  alone  rotary  balance  model,  Eqn  (40)  is  used  to  determine  the 
aerodynamic  coefficients  as  a  function  of  the  body  rates  of  the  aircraft,  angle  of 
attack  and  sideslip  angle.  The  method  of  homotopy  was  incorporated  to  determine 
a  starting  equilibrium  solution.  A  rudder  sweep  was  performed  with  the  McDonnell 
model  using  the  equilibrium  solution  (a,0,6e)  =  (5.0,12.02,-1.39)  identified  by 
Baumann  (7:34).  Because  of  the  dissimilarities  between  the  two  coefficient 
databases,  the  equilibrium  solution  used  for  the  homotopy  application  needed  to 
be  located  at  high  angle  of  attack,  a  >  60°  and  with  zero  control  deflections.  Using 
the  method  of  homotopy  and  the  known  solution  from  the  McDonnell  model  the 
path  of  the  equilibrium  solution  was  transitioned  to  a  stand  alone  rotary  balance 
model  using 

C,  =  "  (1  -  (‘'1) 

where  i=x,y,z,l,m  and  n.  The  Cj  RB  coefficients  are  representative  of  Eqns  (19)  to 
(24).  The  parameter  X  was  used  as  the  control  parameter  in  AUTO  to  trace  the 
path  of  equilibrium  solutions  until  \=1  where  the  model  was  completely  defined  by 
the  rotary  balance  coefficient  data.  This  model  will  be  referred  to  as  the  Rotary 
Balance  (RB)  Model. 

Because  of  the  questionable  quality  of  the  axial  and  side  force  rotary  balance 
coefficients,  a  second  stand  alone  rotary  balance  model  was  created  which  used 
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the  Baumann  model  and  Cy  coefficients  in  lieu  of  the  associated  RB 

coefficients.  Again  using  the  method  of  homotOpy  for  the  remaining  four 
coefficients  a  starting  equilibrium  solution  was  determined.  This  model  will  be 
referred  to  as  the  Modified  Rotary  Balance  (RB)  Model.  The  Modified  RB  Model 
will  assist  in  the  determination  of  the  impact  to  the  investigation  of  the  undesirable 
coefficients.  The  validity  of  combining  the  two  databases  may  be  of  some  concern 
due  to  the  different  testing  conditions  however  results  should  provide  information 
on  the  impact  of  a  perturbed  coefficient  database  and  specifically  the  contributions 
of  the  axial  and  side  force  contribution  to  the  equilibrium  solution  paths. 


When  combining  the  two  coefficient  databases  it  is  important  not  to  duplicate 
information.  The  rotary  aerodynamic  data  contains  static  information.  The  forced 
oscillation  coefficient  contributions  of  the  1988  Aerobase  contain  dynamic 
information  which  is  at  least  partially  duplicated  by  rotary  data.  To  eliminate  the 
duplicate  information,  the  rotary  balance  data  will  be  biased  to  eliminate  static 
contribution.  The  rate  components  p(,5^.,q5,s<.,and  are  used  to  isolate  the  and 
eliminate  the  duplicate  rotary  contributions  of  the  1988  Aerobase.  The  new 
coefficients  are  now  represented  by 

^i.Hybrid  ~  ^i.RB  ^i.RB  Static  ^  ^i.Mc  non-rot 


(42) 


where  i=x,y,z,I,m  and  n.  The  static  rotary  aerodynamic  coefficients  are  determined 
for  a  given  configuration  by  setting  f2ss  =  0-  The  static  and  dynamic  terms  from  the 
1988  Aerobase  coefficients  with  the  rotary  contribution  eliminated  are  represented 
by 

^i,Mc  non-tol  ~ 


Since  the  integrity  of  static  configurations  of  the  rotary  aerodynamic  data  is 
questionable,  a  Modified  Hybrid  model  has  also  been  developed.  The  Modified 
Hybrid  Model  replaces  the  static  contribution  of  the  1988  Aerobase  \A/ith  the  static 
contribution  inherent  in  the  rotary  balance  data.  This  will  investigate  the  strength 
and  weakness  of  the  static  aspect  of  the  rotary  balance  data.  The  Modified  Hybrid 
Model  coefficients  are  defined  by  the  relation 

^/.Modified  Hybrid  ~  ^i,Mc  non-rot  ^i.Mc  Static  ^  ^i.RB 


where  i=x,y,z,l,m,  and  n.  The  static  contributions  of  the  1988  Aerobase,  CjMcsta-ic. 
are  determined  for  a  given  configuration  by  setting  the  body  rates  equal  to  zero. 

Once  the  databases  are  combined,  the  transition  from  the  stand  alone  1988 
Aerobase  and  the  Hybrid  Models  across  a  =  30"  needed  to  be  smoothed.  Figure 
5.4  of  chapter  V  show  the  large  difference  that  could  occur  in  coefficient  values. 
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The  high  and  low  a  databases  are  blended  together  with  the  relation 


Where  Blend  = 


Q  =  ^  -  C;.  J(3  -  2Blend){Blend)^ 

a-30  (45) 


for  angles  of  attack  from  30*  to  35*.  Due  to  the  limitation  of  the  rotary  balance 
database  the  comparison  of  the  models  for  the  analysis  will  be  limited  to  angles 
of  attack  above  30*.  Unlike  previous  studies  mentioned  in  Chapter  1  utilizing  RB 
data  in  bifurcation  analysis  models,  the  span  of  angle  of  attack  was  not  restricted 
by  spin  tunnel  testing  limitations  of  a  >  55°.  Eqn  (45)  was  utilized  as  a  method  of 
homotopy  for  transition  from  known  McDonnell  Model  solutions  in  low  a  to  the 
Hybrid  Model  equilibrium  solutions  in  high  a. 
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V.  Results 


The  objective  of  this  investigation  was  to  deternnine  the  effectiveness  of  using 
rotary  balance  data  for  modeling  of  high  AOA  aircraft  dynamics.  The  best 
indication  of  the  capability  of  a  model  is  to  compare  its  output  to  experimental  flight 
test  data.  For  asymmetric  configurations,  flight  tests  of  the  F-15  showed  highly 
oscillatory  spinning  motion  at  a=50°  to  65°  with  yaw  rates  of  40  to  90  7sec.  These 
results  can  v?7  as  much  as  ±  20°  in  AOA  and  ±  20  in  degrees/sec  in  yaw  rate  for 
a  given  configuration.  Smooth  spin  modes  with  symmetric  loadings  exhibited 
average  AOAs  from  65°  to  75°  with  average  yaw  rates  of  75  to  1 33°/sec.  Spins 
of  higher  rates  are  possible.  However,  because  of  the  pilot’s  physical  limitations, 
flight  testing  does  not  pursue  determination  of  the  maximum  spin  rate  capability 
(3:27-28).  A  typical  full-scale  aircraft  flight  test  showed  a  right  spin  of  a  =  75°  with 
a  spin  rate  of  3  sec/turn  (4:12).  During  spin  tunnel  testing,  to  investigate  free  spin 
modes,  the  scale  model  was  allowed  to  rotate  freely  when  subjected  to  the  free 
stream  air.  The  tests  with  pro-spin  controls,  recovery  controls  and  symmetric 
stabilator  deflection  gave  results  comparable  to  the  full  scale  flight  test  results 
(4:12).  To  investigate  the  effectiveness  of  RB  data,  bifurcation  diagrams 
comparable  to  full  scale  flight  test  results  were  developed  and  then  used  to 
compare  the  aircraft  models  developed. 

Utilizing  equilibrium  solutions  acquired  from  the  McDonnell  model  and  methods 
of  homotopy,  control  surface  variations  were  made  on  each  of  the  five  models:  RB, 
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Modified  RB,  Modified  Hybrid,  Hybrid  and  the  McDonnell.  AUTO  was  used  to 
acquire  equilibrium  starting  points  for  high  a  rudder  deflections,  elevator  deflections 
and  aileron  deflections  representative  of  pro-spin  and  recovery  control 
configurations.  To  acquire  specific  configurations,  alternating  variations  of  different 
control  surfaces  occasionally  had  to  be  exercised  to  acquire  equilibrium  solutions 
on  a  desired  solution  branches  in  the  high  a  regime.  All  bifurcation  diagrams  were 
developed  with  a  fixed  thrust  level  of  8300  Ibf.  The  results  of  the  investigation  will 
be  discussed  in  the  order  of  the  'control  surface  deflections  indicated  with 
numerical  and  graphical  results  presented  at  the  conclusion  of  the  chapter.  The 
bifurcation  diagrams  presented  do  not  always  contain  all  possible  equilibrium 
branches  due  to  the  difficulty  in  obtaining  starting  equilibria  for  unknown  solution 
branches. 

Rudder  Sweep 

As  an  initial  comparison  of  the  rotary  balance  and  Hybrid  models  a  rudder 
sweep  was  performed  at  high  AOA  to  determine  how  different  the  equilibrium 
solution  paths  were.  During  full  scale  flight  tests  it  was  found  that  a  rudder  roll 
entry  technique  was  the  easiest  method  to  intentionally  enter  a  spin.  A  rudder  roll 
would  be  initiated  after  the  aircraft  was  positioned  above  20°  AOA  (3:26).  Similarly 
it  was  identified  by  Baumann  (6:47)  and  McDonnell  (24:30),  an  effective  method 
for  obtaining  an  initial  high  AOA  attitude  was  through  an  elevator  deflection 


64 


followed  by  a  rudder  sweep  (i.e.  continuation)  to  obtain  high  AOA  spin  conditions. 
To  further  investigate  the  development  of  the  equilibrium  solutions,  rudder 
continuations  were  made  at  elevator  deflections  of  -5“,  -19°  and  -25°.  The 
equilibrium  solution  paths  of  the  rudder  continuations  are  shown  in  Figure  5.1  for 
the  five  models. 

Comparison  of  the  RB  Model  and  Modified  RB  Model  (Figure  5.1  a,b)  show  minor 
changes  due  to  variation  in  the  axial  and  side  force  coefficients.  The  RB  model 
is  exhibiting  an  expanding  center  phenomenon  in  the  range  of  a  =  30°  to  40°. 
Cause  of  such  behavior  has  not  been  determined  however  it  may  be  reflective  of 
the  polynomial  behavior  inherent  in  the  RB  coefficient  model. 

The  Hybrid  model  (Figure  5.1c)  shows  similarity  to  the  McDonnell  ModeL(Figure 
5.1  e)  in  the  upper  branch  for  a  ~  85°  with  rudder  deflections  from  -20°  to  10°. 
Otherwise  the  results  are  rather  different.  The  Modified  Hybrid  Model  (Figure  5.1  d) 
shows  indication  that  the  static  contributions  of  the  1988  Aerobase  are  a  major 
contributor  to  the  character  of  the  equilibrium  solutions.  The  Modified  Hybrid 
Model  with  RB  static  contributions  shows  results  closer  to  the  isola  of  the  RB 
Model  as  well  as  the  upper  branch.  Examining  the  development  of  the  solution 
paths  of  the  Hybrid  model  shows  a  closed  path  for  5e=-5°  indicating  possible 
formation  of  an  isola  center.  This  is  indicative  of  the  RB  model  however  the  event 
is  occurring  across  high  a. 

Full  scale  flight  test  results  showed  rudder  deflection,  either  pro-spin  or  anti-spin 
had  no  apparent  effect  on  spin  recovery  (3:28).  Comparing  the  behavior  of  the 
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solution  branches  to  full  scale  flight  tests  shows  the  RB  model  (Figure  5.1a)  having 
the  best  correlation.  The  ineffectiveness  of  rudder  in  a  pro-spin  configuration  at 
high  AOA  with  elevator  stick  aft  is  reflected  in  the  RB  Model  results  with  a 
continuous  branch  for  full  ±  5,  with  no  bifurcation  branches  to  lower  AOA.  The  RB 
model  indicates  large  jumps  in  state  would  need  to  occur  for  transition  to  the 
solution  branches  in  lower  a.  The  McDonnell  model  (Figure  5.1  e)  exhibits  effective 
rudder  deflection  spin  recovery  solution  paths  from  high  to  low  a.  Small  jumps  in 
states  could  occur  as  the  aircraft  states  transition  between  continuation  branches. 
McDonnell  investigated  the  effectiveness  of  rudder  in  his  model  for  spin  recovery. 
He  found  the  rudder  effectiveness  was  caused  by  the  curve  fitting  of  the  stability 
derivative  resulting  in  twice  the  effectiveness  of  the  rudder  at  larger  negative 
values  of  5^.  He  made  the  conclusion  that  the  model  was  somewhat  inaccurate 
in  the  effects  of  rudder  at  high  AOA  (29:39).  This  could  explain  the  many 
differences  found  in  Figure  5.1 . 

The  blending  of  the  rotary  balance  data  with  the  McDonnell  model  in  the  Hybrid 
model  may  not  be  indicative  of  the  capabilities  of  the  Hybrid  model  because  of  the 
identified  problem  with  the  McDonnell/Baumann  coefficient  database.  However, 
it  would  be  expected  for  the  Hybrid  model  to  exhibit  behavior  that  lie  between  the 
two  contributing  parts.  Figure  5.1  (c)  is  exhibiting  unique  behavior  especially  along 
the  high  a  right  spin  branch.  There  may  be  a  detached  branch  that  was  not 
acquired  that  accounts  for  the  solutions  above  TO”  AOA. 
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Elevator  Sweep 


The  rudder  deflection  parameter  sweeps  with  fixed  elevator  and  aileron 
deflections  represented  asymmetric  configurations.  A  second  set  of  bifurcation 
diagrams  were  developed  for  comparison  of  a  symmetric  configuration.  Aileron 
and  rudder  were  held  fixed  at  zero  deflection  while  the  elevator  was  varied  from 
-25'’  to  25*.  Even  though  the  RB  coefficient  data  was  not  defined  for  elevator 
deflections  greater  than  zero,  it  was  felt  the  symmetry  of  the  configuration  could 
support  meaningful  results  above  =  0*. 

Again  comparison  of  the  results  from  the  Modified  RB  model  (Figure  5.2b) 
exhibited  minor  variation  in  state  variables  and  equilibrium  solutions  of  the  RB 
model  (Figure  5.2a).  It  was  found  for  each  configuration  analyzed  during  this 
investigation  the  Modified  RB  Model  exhibited  minor  variations  in  state  variables 
and  equilibrium  solution  branches.  The  results  indicated  that  perturbations  in  axial 
and  side  force  coefficients  had  a  minor  impact  on  the  equilibrium  solutions. 

Again  the  results  of  the  RB  model  (Figure  5.2a)  and  McDonnell  model  (Figure 
5.2e)  are  very  different.  They  both  have  a  branch  of  right  flat  spin  modes  across 
high  a  but  the  RB  model  exhibits  very  linear  behavior.  The  RB  model  was  also 
unable  to  identify  stable  solutions.  Mehra  and  Carroll  (26:76)  found  it  was  usually 
not  possible  to  obtain  flat  spin  equilibria  when  only  static  and  forced  oscillation 
data  are  used.  The  Baumann  and  McDonnell  model  were  able  to  Identify  flat  spin 
equilibria  leading  one  to  believe  that  the  blended  flight  test  data  may  have  been 
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a  major  contributor  to  the  character  of  the  equilibrium  solutions.  The  stable 
regions  identified  by  the  McDonnell  model  are  small  however.  The  results  of  the 
RB  model  follow  the  findings  of  Mehra  and  Carroll  for  no  stable  equilibria.  The 
RB  model  did  exhibit  unstable  branches  in  similar  location  of  the  periodic  wing  rock 
behavior  of  the  McDonnell  model  however  at  a  much  higher  a  and  non-periodic. 

Figures  5.2c  and  5.2e  show  the  Hybrid  model  is  very  similar  in  appearance  to 
the  McDonnell  model.  The  Hybrid  model  result  shows  a  more  drastic  hysteresis 
with  jump  phenomena  very  likely  for  positive  and  negative  elevator  deflections  in 
the  range  8°  to  15*.  Considering  the  concept  of  expanding  centers,  the  Hybrid 
solution  appears  to  be  at  a  different  stage  of  development  compared  to  the 
McDonnell  model.  Again  another  parameter  influence  appears  evident  to  reflect 
the  stages  of  isola  development  between  models.  Between  8^  =  ±  5*  the  a  -  8^ 
bifurcation  diagrams  appear  approximately  the  same.  Table  IV  compares  the 
asr  raft  states  for  8g  =  0*  to  determine  how  different  the  models  are.  The  Hybrid 
Model  is  indicating  slower  rates  resulting  in  a  much  slower  spin.  To  investigate 
the  possible  cause  of  the  differences  of  the  models,  the  force  and  moment 
coefficients  and  states  of  the  aircraft  in  the  region  -25°  <  8^  <  -5°  on  the  high  a 
branch  of  Figure  5.2  were  examined. 

Figure  5.3  shows  the  comparison  of  each  of  the  rernaining  seven  states.  In 
each  diagram  the  states  have  similar  behavior  in  each  model  until  the  turning  point 
at  8fe  =  -15°;  i.e.  the  Hybrid  and  McDonnell  branches  have  similar  states  for 
8e>-1 5°,  until  the  turning  point  on  the  Hybrid  model  at  8e  =  -1 5°  after  which  the 


Table  IV  Comparison  of  right  spin  states  of  the  Hybrid  and 
McDonnell  model  at  high  AOA  with 
Sa=0  • ,  5d=0  • ,  5<,=0  •  and  6^=0*. 


Hybrid  Model 

McDonnell  Model 

a 

66.35* 

64.29* 

P 

-11.48* 

-6.71* 

P 

0 . 4728  rad/sec 

0.7629  rad/sec 

q 

-0.1261  rad/sec 

-0.0901  rad/sec 

r 

1.072  rad/sec 

1.593  rad/sec 

0 

-23.65* 

-25.56* 

(j) 

-6.71* 

-3.236* 

Vtr 

238  ft/sec 

238  ft/sec 

Spin  Rate 

6.36  sec/turn 

3.56  sec/turn 

model’s  exhibit  quite  different  behavior. 

As  discussed  in  Chapter  III,  a  jump  event  could  occur  within  the  hysteresis.  The 
jumps  in  states  would  be  evident  of  abrupt  changes  in  the  aircraft  attitude  and 
rates,  if  the  Hybrid  model  were  a. prediction  of  actual  F-15B  flight  behavior,  the 
aircraft  would  experience  sudden  changes  in  attitude  if  the  elevator  were  held 
between  a  -8°  and  -IS"  elevator  deflection  during  a  spin. 

The  examination  of  the  states  has  given  additional  views  of  the  differences 
between  the  models  however  the  coefficients  are  now  examined  to  assist  in 
identifying  the  cause  of  the  difference.  Figure  5.4  are  plots  of  the  three  force  and 
three  moment  coefficients  values  as  the  elevator  is  swept  between  -5°  and  -25°. 
It  should  be  noted  that  the  coefficients  for  each  model  were  determined  following 
the  locus  of  equilibrium  states  from  Figure  5.3.  The  hysteresis  event  is  very 
evident  in  each  of  the  coefficients.  With  such  strong  differences  in  coefficient 
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values  it  would  be  expected  to  have  very  different  equilibrium  solutions.  The 
bifurcation  diagrams  of  Figures  5.2(c)  and  (e)  may  have  visual  similarity,  but  the 
underlying  state  of  the  aircraft  dynamics  does  not. 

It  would  be  best  to  compare  the  coefficients  at  a  single  set  of  state  values.  The 
force  and  moment  coefficients  were  determined  for  a  fixed  aircraft  state  with  varied 
elevator  deflections  of  -5°  to  -25°.  The  fixed  state  chosen  was  the  flat  right  spin 
exhibited  by  the  McDonnell  model  along  the  upper  branch  of  Figure  5.2e  at  an 
elevator  deflection  of  -15°.  The  aircraft  state  is  exhibited  in  Table  V. 

Table  V  Representative  aircraft  state  for  comparison  of 
the  Hybrid  and  McDonnell  model  coefficient 
behavior . 

a  p  p  q  r  e  (j)  v„ 

68  -5.44  0.69  -0.079  1.714  -21.92  -2.65  236.5 


It  should  be  noted  that  for  bifurcation  analysis  the  diagrams  are  traces  of  the 
equilibrium  solutions  of  the  system.  The  state  in  Table  V  from  the  McDonnell 
model  is  only  an  equilibrium  solution  for  that  model  and  does  not  represent  an 
equilibrium  solution  for  the  Hybrid  model.  Therefore  it  should  not  be  compared 
relative  to  the  behavior  exhibited  in  Figures  5.2  and  5.3. 

Figure  5.5  shows  the  comparison  of  the  coefficients  for  the  aircraft  state  defined 
in  Table  V.  Except  for  the  axial  force  coefficient  (Figure  5.5a) ,  which  is  practically 
the  same,  alLthe  coefficients  appear  to  be  offset  by  a  "static"  bias.  There  are 
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slight  changes  in  slope  along  some  of  the  coefficients  however  the  bias  is  the  most 
predominate  effect.  It  is  most  probable  that  the  questioned  static  aspects  of  the 
modeled  rotary  balance  data  has  introduced  an  error  bias  in  the  Hybrid 
coefficients.  The  Modified  Hybrid  model  (Figure  5.2d)  shows  very  different 
behavior  when  the  static  contributions  are  representative  of  only  the  rotary  balance 
data.  The  dramatic  change  in  the  equilibrium  surface  is  most  likely  caused  by  the 
erroneous  static  data.  If  there  are  errors  in  the  static  data,  as  it  is  either  added  or 
subtracted  from  a  model  it  may  be  removing  essential  information  leaving  noise 
vice  relevant  information. 

Table  VI  is  a  comparison  of  each  of  the  five  models  to  the  spin  tunnel  spin 
mode  predictions  and  the  full-scale  flight  test  spin  mode  of  75“  AOA  and  3  sec/spin 
rate  all  at  common  control  settings.  Comparing  the  Hybrid  model  to  the  spin 
tunnel  predictions,  the  inclusion  of  RB  data  has  not  modeled  the  aircraft  states 
better  than  the  McDonnell  model.  The  Hybrid  model  results  did  not  lie  between 
the  RB  and  McDonnell  model  as  would  be  expected.  There  appears  to  be  an 
additional  factor  influencing  its  results.  The  RB,  Modified  RB  and  Modified  Hybrid 
models  showed  very  good  correlation  to  full  scale  flight  test  results  and  the  spin 
tunnel  predictions. 

Aileron  Sweep 

Various  sweeps  of  different  control  surfaces  were  performed  with  AUTO  to  set 
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up  the  conditions  for  the  spin  modes  identified  in  Tables  V!l  and  VIII.  Once  a 
solution  for  a  specific  configuration  was  acquired,  an  aileron  continuation  was 
made.  By  fixing  all  controls  except  aileron,  the  aircraft  behavior  was  examined  as 
aileron  was  deflected  through  the  desired  pro-spin  or  spin  recovery  configuration. 
Only  right  spins  were  identified. 

Pro-spin  controls  are  acquired  by  obtaining  full  cross-controlled  lateral- 
directional  inputs  (i.efull  positive  aileron  and  full  negative  rudder)  during  an-abrupt 
full  negative  elevator  deflection  (full  aft  stick  position).  Immediately  after  high  a 
attitude  is  obtained,  the  elevator  is  returned  to  a  neutral  position  as  the  aircraft 
settles  into  the  spin  (3:27).  Analogous,  using  AUTO,  alternating  rudder  and 
elevator  sweeps  were  made  to  acquire  high  a  solutions  and  then  an  aileron  sweep 
was  made  to  acquire  a  cross-control  configuration.  A  final  elevator  sweep  was 
made  to  obtain  a  neutral  elevator  deflection.  Using  the  resulting  equilibrium 
solution  an  aileron  continuation  was  performed  resulting  in  Figure  5.6. 

Full  scale  flight  tests  showed  satisfactory  spin  recovery  can  be  accomplished 
with  near  full  aileron/differential  stabilator  deflection  into  the  spin  direction  (3:28). 
For  recovery  from  the  right  spin  identified  in  Figure  5.6;  full  negative 
aileron/differential  elevator  deflection  was  made.  Since  rudder  deflection  has 
already  been  determined  earlier  to  be  ineffective  for  spin  recovery,  its  position  is 
not  critical.  Figure  5.6  also  showed  recovery  control  effectiveness  as  the  aileron 
was  deflected  into  the  right  spin  direction  (negative  deflection).  Figure  5.7  was 
similarly  developed  for  comparison  of  the  specific  data  in  Table  Vlli. 
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The  RB  model  showed  good  correlation  to  the  spin  tunnel  predictions  and  full 
scale  flight  test  predictions  for  the  behavior  in  a  pro-spin  and  recovery  control 
configuration.  The  conditions  identified  in  Table  VII  are  for  pro-spin  Oontrols.  The 
RB  model  (Figure  5.6a)  identifies  the  equilibrium  solution  along  a  near  horizontal 
branch  showing  no  immediate  spin  recovery  to  lower  a  with  moderate 
perturbations  about  83  =  20°.  The  conditions  in  Table  VIII  are  for  recovery  controls. 
The  RB  model  (Figure  5.7a)  identifies  the  equilibrium  solution  along  a  vertical 
branch  indicating  spin  recovery  to  lower  a,  i.e.  recovery  will  occur  as  83  decreases 
towards  -20°. 

Comparing  the  RB  model  between  Figures  5.6  and  5.7  does  provide  evidence 
of  some  effects  of  rudder  on  recovery  controls.  Figure  5.6  with  negative  rudder 
deflection  shows  recovery  with  less  aileron  deflection  positioned  into  the  spin 
direction  than  with  full  positive  rudder  deflection  as  shown  in  Figure  5.7. 

In  contrast  to  the  RB  model  results,  the  McDonnell  model  (Figures  5.6e  and 
5.7e)  did  riot  exhibit  pro-spin  nor  recovery  solutions  for  the  same  configurations. 
McDonnell  noted  the  ineffectiveness  of  aileron  deflection  for  spin  recovery  in  his 
model  however  did  not  pursue  investigating  its  cause  other  than  noting  it  was  a 
cause  of  the  spin  characteristics  of  the  aircraft  and  not  a  problem  v;ith  the  model 
(24:40). 
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Discussion 


The  states  extracted;  for  a  specific  aircraft  configuration  presented  in  Tables 
VI, VII  and  VIII  may  be  deceiving  as  to  the  similarity  between  models.  The  discrete 
examples  present  similar  results  however  the  qualitative  structure  of  the  bifurcation 
diagrams  indicate  very  different  behavior. 

Considering  that  nearly  all  solutions  obtained  were  unstable,  it  was  difficult  to 
justify  that  a  particular  behavior  may  occur.  Unfortunately,  with  the  gross 
differences  between  models,  it  was  hard  to  draw  any  conclusions  to  the  "influence" 
of  rotary  balance  data.  Conclusions  could  be  drawn  however  to  what  information 
it  could  and  could  not  provide  when  examined  by  itself.  It  was  identified  by  the 
results  of  the  Hybrid  and  Modified  Hybrid  model  that  the  static  rotary  balance  data 
contributions  had  a  major  influence  on  the  proper  blending  rotary  aerodynamic 
information  into  the  Hybrid  model.  If  large  errors  were  evident  in  the  static  data 
they  could  cause  meaningful  rotary  information  to  be  "erroneously  removed" 
leaving  noise  perturbations.  This  may  account  for  the  seemingly  independent 
character  of  the  Hybrid  model  with  the  exception  of  elevator  deflection  bifurcation 
diagrams.  The  rotary  balance  data  apparently  has  less  sensitivity  to  the  elevator 
control  surface  than  the  McDonnell  model.  The  model  conflicts  that  occurred  as 
rudder  arid  aileron  deflections  were  made  provided  additional  difficulty  for  the 
Hybrid  model  to  provide  meaningful  information. 
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Table  VI  Clockwise 

spin  state. 

5,=0*  ,5 

,=0-,8e=-25' 

/ 8r=0 ' 

a 

(deg) 

Spin  Rate  f^sb/2Vcr 
(sec/turn ) 

( rad/sec ) ( ft/sec ) 

RB  Model 

75.3 

3.2 

0.16 

1.973 

260 

Modified 

RB  Model 

75.4 

3.2 

0.16 

1.972 

261 

Hybrid 

84.9 

1.7 

0.30 

3.657 

257 

Modified 

Hybrid 

76.5 

5.1 

0.11 

1.220 

245 

McDonnell/ 

Baumann 

82.2 

2.1 

0.29 

3.013 

226 

Spin  Tunnel 
Predict 

65.0 

5.1 

0.10 

1.271 

272 

Flight  Tests 

75.0 

3.0 

-  -  -  - 

2.0 

— 

Table  VII  Pro-spin 

8a=.20* 

controls  model  comparisons. 

=6*  ,5e=0',5r=-15' 

a 

(deg) 

Spin  Rate  £^sb/2Vcr 
(sec/turn) 

as  v„ 

(rad/sec) (ft/sec) 

RB  Model 

84.4 

1.98 

0.25 

3.177 

263 

Modified 

RB  Model 

83.6 

2.04 

0.26 

3.080 

255 

Hybrid 

63.2 

5.28 

0.11 

1.189 

236 

Modified 

Hybrid 

85.9 

1.31 

0.34 

4.794 

301 

McDonnell/ 

Baumann 

60.3 

3.79 

0.15 

1.658 

241 

Spin  Tunnel 
Predict 

80.0 

2.70 

0.21 

2.390 

244 

Flight  Tests 

70-85 

2-6 

— 

1-3 

— 

75 


Table  VIII  Recovery  controls  model  comparisons. 
5^=-20‘  ,5<j=-6* /6e=0‘,5,.=30* 


a 

(deg) 

Spin  Rate  Dssb/2Vcj. 
(sec/turn) 

Pss 

(rad/sec) 

v„ 

(ft/sec) 

RB  Model 

62.2 

4.60 

0.12 

1.363 

250 

Modified 

RB  Model 

62.2 

4.78 

0.11 

1 . 313 

257 

Hybrid 

66.0 

5.17 

0.11 

1.2158 

244 

Modified 

Hybrid 

52.1 

7 .02 

0.08 

0.895 

252 

McDonnell/ 

Baumann 

64.8 

3.50 

0.16 

1.790 

240 

Spin  Tunnel 
Predict 

56.0 

6.10 

0.08 

1.039 

278 

Flight  Tests 

50-65 

4-6 

_  .  _  . 

0.7-1. 6 

— 

(a)  Rotary  Balance  Model. 


(b)  Modified  Rotary  Balance  Model. 
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Figure  5.1  Rudder  Sweep  Bifurcation  Diagram  (5a=0  ,5^=0  ,5e=-5  ,-19  ,-25  ) 
(c)  Hybrid  Model. 


Figure  5.1  Rudder  Sweep  Bifurcation  Diagram  (5a=0°,5j.=0',5e=-5*,-l9°,-25  ) 
(d)  Modified  Hybrid  Model. 
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Figure  5,1  Rudder  Sweep  Bifurcation  Diagram  (5a=0°,5c=0°,5g=-5',-l9*,-25®) 
(e),  McDonnell  Model. 
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Figure  5.2  Elevator  Sweep  Bifurcation  Diagram  (5a=0“,5j;=0“,5,=0°) 
(a)  Rotary  Balance  Model. 
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Figure  5.2  Elevator  Sweep  Bifurcation  Diagram  (5a=0°,5c=0°,5,=0'’) 
(b)  Modified  Rotary  Balance  Model. 
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Figure  5.2  Elevator  Sweep  Bifurcation  Diagram  (5a=0  ,5^=0  ,5^=0  ) 
(c)  Hybrid  Balance  Model. 


Figure  5.2  Elevator  Sweep  Bifurcation  Diagram  (83=0  ,5^=0  ,5,=0  ) 
(d)  Modified  Hybrid  Model. 
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Figure  5.3  Comparison  of  Hybrid  and  McDonnell  Model  Aircraft  States  as  Sg 
is  Varied  (8g=0“,5d=0“,5,=0“).  (a)  Sideslip,  p. 


Figure  5.3  Comparison  of  Hybrid  and  McDonnell  Model  Aircraft  States  as  6^ 
is  Varied  (8g=0“, 5^=0“ ,5=0*).  (b)  Roll  Rate,  p. 


83 


Figure  5.3  Comparison  of  Hybrid  and  McDonnell  Model  Aircraft  States  as  6e 
is  Varied  (5a=0°,5d=0°,5,=0“).  (c)  Pitch  Rate,  q. 


Figure  5.3  Comparison  of  Hybrid  and  McDonnell  Model  Aircraft  States  as  % 
is  Varied  (5a=0°,5^)=0*,5f=0°).  (d)  Yaw  Rate,  r. 
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Figure  5.3  Comparison  of  Hybrid  and  McDonnell  Model  Aircraft  States  as  Sg 
is  Varied  (5a=0°,5d=0°,5,=0'').  (e)  Bank  Angle,  (J). 


Figure  5.3  Comparison  of  Hybrid  and  McDonnell  Model  Aircraft  States  as  Sg 
is  Varied  (53=0°, 5^=0°, 5=0’).  (f)  Pitch  Angle,  0. 
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Figure  5.3  Comparisdh  of  Hybrid  and  McDonnell  Model  Aircraft  States  as  5^ 
is  Varied  (6a=0°.5d=0',6,=0°).  (g)  True  Velocity,  V„. 
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Figure  5.4  Comparison  of  Aerodynamic  Coefficients  as  5^  is  Varied 
(53=0^5d=0^5,=0°).  (a)  Axial  Force  C,. 


Figure  5.4  Comparison  of  Aerodynamic  Coefficients  as  5^  is  Varied 
(53=0^5d=0^6  =0").  (b)  Side  Force  Cy. 


Figure  5.4  Comparison  of  Aerodynamic  Coefficients  as  5e  is  Varied 
(S3=0^6d=0^8,=0“).  (c)  Normal  Force  C,. 
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Figure  5.4  Comparison  of  Aerodynamic  Coefficients  as  5e  is  Varied 
(83=0“, 6^=0°, 5  =0“).  (d)  Roil  Moment  C,. 
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Figure  5.4  Comparison  of  Aerodynamic  Coefficients  as  5^  is  Varied 
(5^=0°, 5^=0*. 5, =0’’).  (e)  Pitch  Moment  C^. 


Figure  5.4  Comparison^of  Aerodynamic  Coefficients  as  5^  is  Varied 
(5a=0“,5^=0°,5,=0’).  (f)  Yaw  Moment  C„. 
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Side  Force  Coefficient 


Elevator  Deflection  (deg) 

Figure  5.5  Compdiison  of  the  Aerofjynamic  Coefficients  as  5e  is  Varied  for  the 
state  defined  in  Table  V.  (a)  Axial  Force 


Figure  5.5  Comparison  of  the  Aerodynamic  Coefficients  as  5^  is  Varied  for  the 
state  defined  in  Table  V.  (b)  Side  Force  Cy. 


Figure  5.5  Comparison  of  the  Aerodynamic  Coefficients  as  is  Varied  for  the 
state  defined  in  Table  V.  (c)  Normal  Force  C^. 
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Figure  5.5  Comparison  of  the  Aerodynamic  Coefficients  as  5e  is  Varied  for  the 
state  defined  in  Table  V.  (d)  Roll  Moment  C,. 


Figure  5.5  Comparison  of  the  Aerodynamic  Coefficients  as  5^  is  Varied  for  the 
state  defined  in  Table  V,  (e)  Pitch  Moment  C^,. 


Figure  5.5  Comparison  of  the  Aerodynamic  Coefficients  as  6^  is  Varied  for  the 
state  defined  in  Table  V.  (f)  Yaw  Moment  Cn. 


-30  -10  -10  0  10  20  30 

Aileron  Defleciion  (deg) 


Figure  5.6  Aileron  Sweep  Bifurcation  Diagram  With  Pro-spin  Controls 
(5<.=0.352,5g=0“,5,=-l5*’)  (a)  Rotany  Balance  Model. 


Figure  5.6  Aileron  Sweep  Bifurcation  Diagram  With  Pro-spin  Controls 

(5o=0.35a,5g=0’,5,=-1 5’)  (b)  Modified  Rotary  Balance  Model. 
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Figure  5.6  Aileron  Sweep  Bifurcation  Diagram  With  Pro-spin  Controls 
(5^=0.383, 5e=0°, 5, =-15”)  (c)  Hybrid  Model. 


Figure  5.6  Aileron  Sweep  Bifurcation  Diagram  With  Pro-spin  Controls 
(5d=0.35a,53=0’,5=-15“)  (d)  Modified  Hybrid  Model. 


Figure  5.6  Aileron  Sweep  Bifurcation  Diagram  With  Pro-spin  Controls 
(5d=0.383,5e=0°.5=-15“)  (e)  McDonnell  Model. 
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Figure  5.7  Aileron  Sweep  Bifurcation  Diagram  With  Recovery  Controls 
(5d=0.35a,5e=0°,5,=30°)  (a)  Rotary  Balance  Model. 


Figure  5.7  Aileron  Sweep  Bifurcation  Diagram  With  Recovery  Controls 
(8^=0.383, 5e=0°, 5, =30°)  (b)  Modified  Balance  Model. 
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Figure  5.7  Aileron  Sweep  Bifurcation  Diagram  With  Recovery  Controls 
(6^=0.383, 8,=0°, 5  =30“)  (c)  Hybrid  Model. 


Aileron  Deflection  fdeg) 

Figure  5.7  Aileron  Sweep  Bifurcation  Diagram  With  Recovery  Controls 
(8d=0.38a,5e=0“,5,=30“)  (d)  Modified  Hybrid  Model. 
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VI.  Conclusions 


As  evidenced  from  this  investigation,  one  of  the  distinctive  differences  and 
sources  of  difficulty  in  models  of  aircraft  dynamics  is  the  representation  of  the 
aerodynamic  forces  and  moments.  It  is  difficult  to  accept  any  research  results  that 
apply  experimental  data  for  modeling  without  a  validation  of  the  integrity  of  the 
data  upon  which  the  research  was  based.  This  investigation  presented  an 
opportunity  to  compare  three  models  based  on  three  very  different  aerodynamic 
coefficient  databases.  The  results  have  given  evidence  of  the  caution  that  needs 
to  be  exercised  when  models  are  compared  as  v;ell  as  the  possible  error  that  is 
introduced  when  two  different  sets  of  data  are  combined. 

Both  the  McDonnell  and  RB  model  have  the  ability  to  predict  high  AOA  behavior 
of  the  F-15B.  The  fundamental  difficulty  is  that  the  qualitative  and  quantitative 
outlay  of  equilibria  is  very  different.  There  is  still  the  question  of  which  model  is 
better.  The  McDonnell  model,  representing  static  and  forced  oscillation  data,  has 
demonstrated  wing  rock  beha''ior  indicative  of  full  scale  flight  test  results  as 
identified  in  reference  (3).  It  does  however  have  problems  with  the  consequences 
of  aileron  and  rudder  deflection  at  high  a  which  is  probably  more  evident  of  the 
modelling  of  the  stability  derivatives  than  with  the  capability  of  static  and  forced 
oscillation  data.  With  the  possible  problems  with  the  McDonnell/Baumann 
coefficient  modelling  in  83  and  5.  it  is  difficult  to  draw  any  conclusions  concerning 
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its  relative  merits  when  compared  to  the  RB  model.  The  McDonnell  model  would 
however  provide  initial  indication  of  high  a  spin  behavior  to  initiate  further 
investigations  with  the  RB  model.  Moreover,  the  McDonnell  model  was  used  at 
6a=5r=0  with  thrust  vectoring  producing  quality  results,  so  it  is  not  aiways  limited 
by  this  problem  (24). 

The  results  of  Chapter  V  have  shown  that  rotary  balance  data  does  enable 
prediction  of  aircraft  spin  motion  with  good  correlation  to  full  scale  flight  test 
results. 

In  particular: 

(1 )  The  RB  model  properly  represented  rudder  motion  as  not  being  optimum  for 
spin  recovery  as  demonstrated  by  full  scale  fligh';  tests. 

(2)  The  RB  model  properly  represented  recovery  control  behavior  with  ailerons 
in  a  representative  flight  test  recovery  configuration. 

(3)  in  pro-spin  control  configuration  representative  of  full  scale  flight  test 
behavior,  the  RB  model  showed  no  immediate  spin  recovery. 

in  addition: 

(4)  The  rotary  balance  model  was  unable  lO  !:Isntify  stoble  equilibrium  branches, 
as  predicted  by  Mehra  and  Carroll  for  fighter  aircraft  and  shown  in  their  research 
results  for  the  F-4  fighter  (26). 

The  results  from  the  RB  model  does  demonstrate  rotary  balaiice  data  has  strength 
in  representing  spinning  motions  of  the  aircraft  however  there  is  not  enough 
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evidence  to  declare  supremacy  over  fhe  McDonne”  model. 

These  results  were  determined  with  a  single  polynomial  representation  of  each 
of  the  aerodynamic  force  and  moment  coe^’  .r !  It  should  be  evident  that  the 
modeling  of  the  RB  data  has  smotihec  n.Vi  the  character  of  the  original 
database.  However,  even  with  a  "smoothed''  '.jpresentation  of  the  RB-;data,  it  was 
able  to  provide  good  correlation  to  experimt*  -light  test  results.  The  promising 
results  shown  by  the  RB  data  as  a  stand  alor.e  i  epreser.tation  of  the  aerodynamic 
coefficients  in  high-  AOA  dynamics  analysis  yields  further  investigation  of  the 
effects  of  the  static  errors.  With  the  RB  data  more  effectively  modelled  it  may 
provide  more  character  in  the  bifurcation  diagrams. 

The  problems  encountered  with  the  blending  of  two  databases  acquired  from 
different  testing  facilities  was  evident  in  the  Hybrid  model.  Because  of  the  different 
sources  of  the  data,  inherent  ^!Tors  may  conflict  or  even  amplify.  It  needs  to  be 
assured  that  the  data  being  blended  is  representative  of  the  same  testing 
condition.  The  error  evident  in  the  RB  coefficient  static  configuration  had  a 
definite  influence  on  the  results  found  with  the  Hybrid  model.  The  problems  with 
aileron  and  rudder  deflection  in  the  McDonneli/Baumann  model  introduced  error 
with  the  Hybrid  model.  Any  research  intending  to  use  a  hybrid  model  needs  to 
develop  three  models.  The  RB,  hybrid  and  static  and  forced  oscillation  model, 
need  to  be  examined  to  ensure  eaci.  is  representing  basic  behavior  such  as  rudder 
ineffectiveness  or  anticipated  recovery  results  before  the  hybrid  is  used.  If  the 
databases  were  blended  without  concern,  the  Hybrid  model  may  be  presenting 
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false  information  which  is  based  on  conflicting  databases. 

Results  from  the  Modified  RB  model  showed: that  perturbation  in  the  Axial  and 
Side  Force  coefficients  had  a  minor  influence  on  the  equilibrium  .ijlution  branches 
when  investigating  the  spin  regime  of  the  F-15B.  The  Modified  Hybrid  model  gave 
indication  that  the  static  contributions  to  the  aerodynamic  coefficients  was  most 
influential  in  the  character  of  equilibrium  solutions. 

Recommendations 

1 .  An  investigation  of  the  evident  static  contribution  problem  of  the  rotary  balance 
data  and  its  effects  on  the  Hybrid  model  could  be  made.  A  filtering  process  could 
be  developed  to  adjust  the  removal  of  the  rotary  balance  static  contribution  at 
small  increments,  artificiaiiy  reducing  its  contribution  to  the  model.  horhotopy 
variable  could  be  defined  to  adjust  the  proportions  of  the  rotary  balance  static 
contributions  removed  and  then  a  continuation  on  this  variable  could  be  performed 
to  see  what  features  of  the  bifurcation  diagrams  are  most  effected.^ 

2.  Ti  a  static  and  forced  oscillation  contributions  of  the  McDonnell  model 
coefficients  could  be  isolated  and  analyzed  using  bifurcation  analysis  to  provide 
insight  into  which  features  of  the  equilibrium  solutions  are  driven  by  static  and 
forced  oscillation  information  separately.  The  results  from  the  analysis  of  the 
McDonnell  model  static  contributions  may  provide  assistance  in  determining  the 
level  of  adjustments  that  could  be  made  on  the  rotary  balance  data  in 


recommendation  (1), 

3.  As  demonstrated  by  the  Modified  RB  model,  analysis  of  the  influence  of  each 
force  and  moment  coefficient  on  the  equilibrium  solutions  should  be  investigated. 
Before  bifurcation  analysis  can  be  used  confidently  for  analysis  of  aircraft  dynamics 
its  tolerance  to  perturbations  in  each  of  the  force  and  moment  coefficients  needs 
to  be  determined.  In  addition,  the  analysis  may  provide  insight  into  the  behavioral 
differences  between  the  models.  Using  the  McDonnell, Rotary  balance  or  an 
artificial  database,  each  coefficient  could  be  perturbed  separately  with  either  a 
static  bias  or  a  selected  function.  ‘l  he  r8s.j.ts  could  assist  in  identification  ot  which 
coefficients  are  driving  the  gross  oifferences.  The  data  could  be  adjusted  to 
develop  certain  phenomena  and  then  by  adjusting-  the  coefficient  data  until  the 
phenomena  distorts  or  disappears,  the  tolerance  in  accuracy  of  the  coefficients 
could  be  determined. 
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Appendix  A:  F-15B  Specification  Data 


The  F-15  specifications  for  the  F-15  are  contained  in  Table  IX.  The  data  was 
obtained  from  Beck  (8)  and  (23). 


Table  IX.  Physical  Characteristics  of  the  F-15B 


Wing 

Area  (theoretical) 

Aspect  Ratio 
Airfoil 
Root 
Xw  155 
Tip 
Span 

Taper  Ratio 

Root  Chord  (Theoretical) 

Tip  Chord 

Mean  Aerodynamic  Chord 

Leading  Edge  Sweep  Angle 

25%  Chord  Sweep  Angle 

Dihedral 

Incidence 

Twist  at  Tip 

Aileron  Area 

Flap  Area 

Speed  Brake  Area 


Control  Surface  Movement 
Aileron 
Speedbrake 
Flap 

Horizontal  Tail 
Rudder 


608  sq  ft 
3.01 

NACA64006.6 

NACA64A(x)04.6  (a  =  0.8  MOD) 
NACA64A203  (  a  =  0.8  MOD) 

42.8  ft 
0.25 

273.3  in 

68.3  in 

191.3  in 
45* 

38.6* 

-1* 

None 

None 

26.5  Sq  ft 

35.8  sq  ft 

31 .5  sq  ft 


±20* 

45*  up 
30°  down 
29*  down,  15*  up 
±30* 
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Vertical  Tail 


Area  (Theoretical  Each) 

62.6  sq  ft 

Rudder  Area  (Each) 

1 0.0  sq  ft 

Span 

10.3  ft 

Root  Chord 

115.0  in 

Tip  Chord 

30.6  in 

Airfoil  -  Root 

NACA0005-64 

-Tip 

NACA0003.5-64 

Taper  Ratio 

0.27 

Leading  Edge  Sweep  Angle 

36.6° 

25%  Chord  Sweep  Angle 

29.7° 

Mean  Aerodynamic  Chord 

81.0  in 

Cant 

2°  out 

Length  (.25c„  to  .25cJ 

212.4  in 

Horizontal  Tail 


Area  (Theoretical) 

120.0  sq  ft 

Area  (Actual) 

1 1 1 .4  sq  ft 

Span 

15.7  ft 

Aspect  Ratio 

2.05 

Taper  Ratio 

0.34 

Root  Chord 

137.2  in 

Tip  Chord 

46.5  in 

Airfoil  -  Root 

NACA0005.5-64 

-Tip 

NACA0002.5-64 

Leading  Edge  sweep  Angle 

50° 

25%  Chord  Sweep  Angle 

43.6° 

Mean  Aerodynamic  chord 

99.3  in 

Dihedral 

0° 

Length  (.25c^^  to  .25Ch) 

241 .0  in 

Wetted  Area 

Fuselage 

1405  sq  ft 

Nozzles 

53  sq  ft 

Horizontal  Tail 

216  sq  ft 

Vertical  tail 

257  sq  ft 

Wing 

698  sq  ft 

Total  Area 

2629  sq  ft 
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Miscellaneous  Data 
Aircraft  Length  63.8  ft 

Aircraft  Height  1 8.6  ft 

Aircraft  Volume  1996  cu  ft 

Aircraft  Gross  Weight  37000  lbs 

C.G.  Station  X  Direction  557.173 

Y  Direction  0.0 

Z  Direction  116.173 


Inertia  Data  is  for  a  basic  F-15  with  4  AIM-7F  missiles,  ammo,  50%  fuel  and  gear 
up. 

lx  25480  slug-fP 

ly  166620  slug-ft^ 

Iz  186930  slug-ft^ 

1x2  -1 000  slug-fP 
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Appendix  B:  Sian  Conventions 


The  airplane  is  considered  in  an  upright  attitude  vWth  all  directions  with  respect 
to  the  pilot  seated  in  the  cockpit.  Refer  to  figure  B-1 .  The  sign  convention  data  for 
the  F-15  was  obtained  from  (23). 


Airplane 

Forces  -  Positive  up,  aft  or  to  the  right 

Moments  -  Positive  when  the  nose  pitches  up,  to  the  right,  or  the  left  wing 
rises. 


Angular 

Velocity  -  Positive  when  the  nose  rotates  up.  to  the  right,  or  the  left  wing  rises. 
Control  Surfaces 

Aileron  -  Positive  when  right  aileron  is  down. 

Differential 

Horizontal  Tail  -  Positive  when  right  panel  is  down. 

Symmetrical 

Horizontal  tail  -  Positive  when  trailing  edge  is  down. 

Rudder  -  Positive  when  trailing  edge  is  to  the  left. 
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Appendix  C:  AUTO  Driver  Program 


C 

C  CAPTAIN  RALPH  D.  FERO  AFIT  GA-91D 
C  MASTERS  THESIS 
C 

C  THE  FOLLOWING  DRIVER  PROGRAM  IS  A  MODIFIED  VERSION  OF 
C  CAPT  ROBERT  MCDONNELL’S  1990  THESIS.  THE  ADDITIONS  IN 
C  THIS  VERSION  OF  THE  PROGRAM  INVOLVE  INTEGRATION  OF 
C  F-1 5  ROTARY  AERODYNAMIC  DATA.  THIS  PROGRAM  SOLVES  THE 
C  NONLINEAR  DIFFERENTIAL  EQUATIONS  OF  MOTION  FOR  THE 
0  F-15B  AIRCRAFT.  FOR  THIS  RESEARCH  EFFORT.  THIS  PROGRAM 

C  WILL  BE  USED  TO  INVESTIGATE  THE  EFFECTS  OF  ROTARY 
0  BALANCE  DATA  ON  THE  ANALYSIS  OF  HIGH  ANGLE  OF  ATTACK 
C  PHENOMENA.  THE  PROGRAM  IS  CAPABLE  OF  VARYING 
0  ELEVATOR.  AILERON.  RUDDER  DEFLECTIONS.  ENGINE 
C  THRUST  VECTOR  (PITCH  AND  YAW).  PORT  AND  STARBOARD 

C  ENGINE  THRUST,  AND  TOTAL  THRUST.  IN  ADDITION.  THE 

C  PROGRAM  HAS  BEEN  MODIFIED  TO  VARY  DIFFERENTIAL 
G  ELEVATOR  AND  A  HOMOTOPY  BLENDING  PARAMETER  FOR 
0  TRANSITION  FROM  THE  1988  F-15  AEROBASE  DATA  TO  ROTARY 
C  AERODYNAMIC  DATA. 

C 

C  LAST  EDITED  ON  14  OCTOBER  1991 
C 

IMPLICIT  DOUBLE  PRECISION(A-H,0-Z) 

COMMON  /RBPOLY/  PC,  OMEGA 
DIMENSION  W(300000),  IW{1000),PC(6,79) 


OPEN(UNIT=3,FILE=’fort.3’) 
OPEN(UNIT=4,FILE=’fort.4’) 
OPEN(UNIT=7,FILE=’fort.7’) 
OPEN(UNIT=8,FILE=’fort.8’) 
OPEN{UNIT=9,FILE=’fort.9’) 
OPEN(UNIT=1 0,FILE=’fort.1  O’) 
OPEN(UNIT=1 2,FILE=’seize’) 
OPEN(UNIT=13.FILE=’seizet’) 

REWIND  7 
REWIND  8 
REWIND  9 
REWIND  10 
REWIND  3 
REWIND  4 
REWIND  12 
REWIND  13 

ADDED  15  AUG  91 
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c 

CALL  RBPOLYCOEF 

CALL  AUTO  -  CONTINUATION  &  BIFURCATION  LOCATION 
SUBROUTINE 

CALL  AUTO(W.IW) 

STOP 
END 


SUBROUTINE  FUNC(NDIM,NPAR,U,ICP.PAR,IJAC,F,DFDU,DFDP) 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  /KS/  K1.K5,K7,K8,K9,K10,K12,K13,K14,K15,K16,K17 
COMMON  /ACDATA/  BWING, OWING, SREF.RHO.RMASS 
DOUBLE  PRECISION  K1,K5,K7,K8,K9,K10,K12,K13,K14,K15,K16,K17 
COMMON  /SEIZE/  CX.CY.CZ,CLM,CMM,CNM 
COMMON  /SEIZET/  CXT.CYT.CZT.CLMT.CMMT.CNMT 
COMMON  /SEIZER/  CXR.CYR.CZR.CLMR.CMMR.CNMR 
COMMON  /RBPOLY/  PC,  OMEGA 


DIMENSION  DFDU(NDIM,NDIM),DFDP(NDIM,NPAR),DELF1  (8), 
+  DELF2(8),U(NDIM),PAR(10),F(NDIM),DX(8),PC(6,79) 


INITIALIZE  SOME  CONSTANTS  THAT  ARE  PASSED  THROUGH 
THE  COMMON  BLOCK  ACDATA 

DATA  IS  FROM  MCAIR  REPORT#  A4172  AND  AFFTC-TR-75-32 
F-15A  APPROACH-TO-STALL/STALL/POST-STALL  EVALUATION 

BWING  -  A/C  WINGSPAN,  FT 

CWING  -  A/C  MEAN  AERODYNAMIC  CHORD,  FT 

SREF  -  A/C  WING  REFERENCE  AREA,  SO  FT 

RHO  -  AIR  DENSITY  AT  20000  FT  ALTITUDE,  SLUG/FT^3 

RMASS  -  A/C  MASS,  SLUGS 

BWING=42.8 

CWING=15.94 

SREF=608. 

RHO=.0012673 

RMASS=37000./32.174 


DETERMINE  CONSTANTS  K1  THROUGH  K17.  SOME  ARE  MADE 
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C  COMMON  AND  PASSED  TO  SUBROUTINE  FUNX  AND  USED  IN  THE 
C  EQUATIONS  OF  MOTION  THERE 
C 

C  INERTIAS  HAVE  UNITS  OF  SLUG-FT'^2 
C 

C  K1  HAS  UNITS  OF  1/FT 
C 

C  K6,  K8.  K11,  K14,  AND  K17  HAVE  UNITS  OF  1/FT''2 

C 

C 

C  IX=  25480.0C10 

C  IY=  166620.0d0 

C  IZ=  186930.0d0 

C  IXZ=  -lOOO.OdO 

C  K1=0.5d0*RHO*SREF/RMASS 
C  K2=(IZ-IY)/IX 

C  K3=IXZ‘IXZ/(IX’IZ) 

C  K4=(IY-IX)/IZ 

C  K5=IXZ/IX 

C  K6=0.5dO*RHO*BWING‘SREF/IX 
C  K7=IXZ/IZ 

C  K8=0;5dO*RHO*SREF*CWING/IY 
C  K9=(IZ-IX)/IY 

C  K10=IXZ/IY 

C  K1 1  =0.5dO*RHO*SREF'BWING/IZ 

C  K12=(K2+K3)/(1.0d0-K3) 

C  K13=(1.0d0-K4)*K5/(1.0d0-K3) 

C  K14=K6/(1.0d0-K3) 

C  K15=(K3-K4)/(1.0d0-K3) 

C  K16=(1.0d0+K2)*K7/(1.0d0-K3) 

C  K17=Kl1/(1.0d0-K3) 

C 

C 

K1  =  3.350088890D-04 
K5  =-3.924646781  D-02 
K7  =-5.3495961 05D-03 
K8  =  3.685650971  D-05 
K9  =.96897131196 
K10  =-6.001 680471 D-03 
K12  =  .79747314581 
K1 3  =-9.61 5755341  D-03 
K14  =  6.472745847  D-04 
K15  =-.754990553922 
K16=  K13 

K17  =  8.822851 558D-05 


FIND  THE  VALUES  OF  F(1)  THROUGH  F(NDIM).  SUBROUTINES 
COEFF  AND  FUNX  ARE  CALLED  ONCE. 
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CALL  COEFF(U,PAR,NDIM,ICP) 

THE  FOLLOWING  ADDED  15  AUG  91 

THE  SUBROUTINE  RBPOLYCOEF  AND  RBCOEF  DETERMINE  THE 
ASSOCIATED  COEFICIENTS  USING  ROTARY  BALANCE  DATA 
OBTAINED  FROM  THE  NASA  LANGLEY  SPIN  TUNNEL  THE  DATA 
WAS  OBTAINED  ON  FLOPPY  DISKS.  THE  DATA  IS  DOCUMENTED 
IN  NASA  CR  3478.  THE  DATA  IS  ONLY  CODED  FOR  AOA  ABOVE 
30  DEGREES  THEREFORE  WILL  BE  BLENDED  AT  THIS  VALUE. 

IF  (U(1)  .GT.  30)  THEN 

CALL  RBCOEF(U,PAR,NDIM) 

SUBROUTINE  BLEND  PERFORMS  THREE  FUNCTIONS.  USING  THE 
PARAMETER  BLEND,  (1)  IDENTIFY  A  UNIQUE  EQUILIBRIUM 
STATE  SOLUTION  FOR  PURE  ROTARY  BALANCE  DATA  BASED 
MODEL, (2)  PERFORM  THE  ADDITION  THRUST  CONTRIBUTIONS  TO 
THE  PURE  ROTARY  BALANCE  DATA  BASED  MODEL,  (3)  PERFORM 
THE  BLEND  TRANSITION  FROM  MCDONNELL’S  MODEL  TO  THE 
HYBRID  ROTARY  BALANCE  DATABASE  MIXED  MODEL  (HYBRID 
MODEL). 

CALLRBBLEND(U,PAR,NDIM,ICP) 

ENDIF 

CALL  FUNX(NDIM,U,F) 


IF(IJAC.EQ.O)  RETURN 
SET  THE  VALUES  OF  DX 

MODIFIED  TO  SCALE  DX  ACCORDING  TO  VARIABLE 
13  JUN  88 


DX0=1.0D-9 

DX(1)=DX0*50.0C10 

DX(2)=DX0*10.0ci0 

DX(3)=DXO*0.5cJ0 

DX(4)=DX0'0.25d0 

DX(5)=DXO*0.5dO 

DX(6)=DX0*50.0d0 

DX(7)=DX0*50.0d0 

DX(8)=DXO‘0.5dO 
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NEXT  THE  PARTIAL  OF  F  W.R.T.  A  GIVEN  PARAMETER  ARE 
FINITE  DIFFERENCED 

PTEMP=PAR(ICP) 

PAR(ICP)=PTEMP+DX(1) 


CALL  COEFF(U,PAR.NDIM,ICP) 


FOLLOWING  ADDED  15  AUG  91 

IF  (U(1)  .GT.  30)  THEN 
CALL  RBCOEF(U,PAR,NDIM) 

CALL  RBBLEND(U,PAR,NDIM,ICP) 
ENDIF 

CALL  FUNX(NDIM,U.DELF1) 
PAR(ICP)=PTEMP-DX(1) 


CALL  COEFF(U,PAR,NDIM,ICP) 


FOLLOWING  ADDED  15  AUG  91 

IF  (U(1)  .GT.  30)  THEN 
CALL  RBC0EF(U,PAR,ND1M) 

CALL  RBBLEND(U,PAR,NDIM,ICP) 
ENDIF 


CALL  FUNX(ND!M,U,DELF2) 

DO  13  l=1,NDIM 

DFDP(l,ICP)=(DELFl(l)-DELF2(l))/(2.0d0*DX(1)) 


13  CONTINUE 
PAR(ICP)=PTEMP 
C 

C  THE  NEXT  DO  LOOP  CALCULATES  THE  PARTIAL  DERIVATIVE  OF 
C  F  W.R.T.  TO  U  USING  FINITE  DIFFERENCES. 

C 

C  SET  U(J)  EQUAL  TO  U+DU,  THEN  CALL  COEFF  WITH  THIS 
C  UPDATED  STATE  VECTOR.  THIS  IS  DONE  SIMILARLY  WITH 
C  U-DU 
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c 

c 

c 


DO  20  J=1,NDIM 
UTEMP=U(J) 
U(J)=UTEMP+DX(J) 


CALL  COEFF(U,PAR,NDIM,ICP) 


FOLLOWING  ADDED  15  AUG  91 

IF  (U(1)  .GT.  30)  THEN 
CALL  RBCOEF(U,PAR,NDIM) 

CALL  RBBLEND(U,PAR,NDIM,ICP) 
ENDIF 


CALL  FUNX(NDIM,U.DELF1) 
U(J)=UTEMP-DX(J) 


CALL  COEFF(U,PAR,NDIM,ICP) 

FOLLOWING  ADDED  15  AUG  91 

IF  (U(1)  .GT.  30)  THEN 
CALL  RBCOEF(U,PAR,NDIM) 

CALL  RBBLEND(U,PAR,NDIM,ICP) 
ENDIF 


CALL  FUNX{NDIM,U,DELF2) 


DO  16  l=1,NDIM 

DFDU(l,J)=(DELF1(l)-DELF2(l))/(2.0d0*DX{J)) 
16  CONTINUE 
C 

U(J)=UTEMP 

C 

20  CONTINUE 

RETURN 

END 
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SUBROUTINE  FUNX(NDIM,U,F) 


SUBROUTINE  FUNX  EVALUATES  THE  NDIM  EQUATIOIm’S  GIVEN  THE 
STATE  VECTOR  U. 

NDIM-  THE  DIMENSION  OF  THE  PROBLEM 
U  -  THE  VECTOR  OF  STATES  ALPHA.  BETA, ...  (INPUT) 

F  -  THE  VECTOR  RESULT  OF  FUNCTION  EVALUATIONS 
(OUTPUT) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

COMMON  /SEIZE/  CX.CY.CZ.CLM.CMM.CNM 
COMMON /SEIZET/CXT,CYT,CZT,CLMT,CMMT,CNMT 
COMMON  /KS/  K1,K5,K7.K8,K9.K10.K12,K13,K14,K15,K16.K17 
DOUBLE  PRECISION  K1,K5,K7,K8,K9,K10.K12,K13.kl4,K15,Kl6,K17 
DIMENSION  U(NDIM),F(ND!M) 

SET  TRIGONOMETRIC  RELATIONSHIPS  OF  THE  STATES  ALPHA, 
BETA,  THETA,  AND  PHI  AND  THEN  SET  P,  Q,  R,  AND  VTRFPS 

IWRITE=1 


DEGRAD=57.29577951  DO 

CA:=COS(U(1)/DEGRAD) 

SA=SIN(U(1)/DEGRAD) 

CB=COS(U(2)/DEGRAD) 

SB=SIN(U(2)/DEGRAD) 

CTHE=C0S(U(6)/DEGRAD) 

STHE=SIN(U(6)/DEGRAD) 

CPHI=COS(U(7)/DEGRAD) 

SPHI=SIN(U(7)/DEGRAD) 

P=U(3) 

Q=U(4) 

R=U(5) 

VTRFPS=1000.0cl0*U(8) 

SET  THE  GRAVITATIONAL  CONSTANT,  FT/SEC 
G=32.1740d0 

THE  FOLLOWING  SYSTEM  OF  NONLINEAR  DIFFERENTIAL 
EQUATIONS  GOVERN  AIRCRAFT  MOTION 

UPDATED  FOR  PROPER  DEGREE-RADIAN  UNITS  AND  PROPERLY 
SCALED  VELOCITY  EQUATION:  7  JUN  88 
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F(1)=ALPHA-D0T 

F(1)=Q+(-(KrVTRFPS*CX-G*STHE/VTRFPS+R'SB)*SA+(Kl*VTRFPS 
+  *CZ+(G*eTHE*CPHI/VTRFPS)-P*SB)'CA)/CB 
F(1)=F(1)'DEGRAD 


F(2)=BETA-D0T 

F(2)=-((K1*VTRFPS*CX-G*STHE/VTRFPS)*SB+R)‘CA+(KrVTRFPS*CY 
+  +G*CTHE*SPHI/VTRFPS)'C8-((K1  ’VTRFPS*CZ+G*CTHE*CPHI/VTRFPS) 
+  *SB-P)*SA 

F(2)=F(2)*DEGRAD 


F(3)=P-D0T 

F(3)=-K12'Q'R4!<13*P'Q+K14-(CLM+K7*CNM)*VTRFPS*VTRFPS 

F(4)=Q-D0T 

F(4)=K8*VTRFPS*VTRFPS*CMM+K9*P'R+K10*(R'R-P*P) 

F{5)=R-D0T 

F(5)=K1 5*P*Q-K1 6*Q'R+K1 7*VTRFPS*VTRFPS'(K5*CLM+CNM) 

F(6)=THETA-D0T 

F(6)=Q*CPHI-R*SPHi 

F(6)=F(6)'DEGRAD 


F(7)=PHI-D0T 

F(7)=P+Q*(STHE/CTHE)*SPHI+R*(STHE/CTHE)*CPHI 

F(7)=F(7)'DEGRAD 


F(8)=VTRFPS-D0T  (SCALED  BY  A  FACTOR  OF  lOOO) 

F(8)=U(8)*((K1  'VTRFPS*CX-G'STHE/VTRFPS)*CA*CB+(K1  ‘VTRFPS'CY 
+  +G*CTHE’SPHI/VTRFPS)*SB 
+  +(K1  ‘VTRFPS*CZ+G*CTHE*CPHI/VTRFPS)'SA*CB) 
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RETURN 

END 

SUBROUTINE  STPNT(NDIM,U,NPAR,ICP,PAR) 


THIS  SUBROUTINE  SETS  THE  VALUES  OF  THE  STATES  AND 
PARAMETERS  AT  THE  START  OF  THE  ANALYSIS.  THE  STATES 
AND  CONTROL  SURFACE  SETTINGS  REPRESENT  AN  EQUILIBRIUM 
STATE  OF  THE  AIRCRAFT 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

C 

DIMENSION  U(NDIM),PAR(10) 

C  U(1)- ALPHA.  DEG 
C  U(2)  -  BETA,  DEG 

C  U(3)  -  P.  RAD/SEC 

C  U(4)  -  Q,  RAD/SEC 

C  U(5)  -  R,  RAD/SEC 

C  U(6)  -  THETA,  DEG 

C  U(7)  -  PHI,  DEG 

C  U(8)  -  TRUE  VELOCITY,  IN  THOUSANDS  OF  FT/SEC 
C 

C  THE  STARTING  POINT  (VECTOR) 

C 

OPEN(UNIT=15,FILE=’fort.15’) 

REWIND  (15) 

C 

READ(15.*)  U(1) 

READ(15,*)  U(2) 

READ(15,’)  U(3) 

READ(15,*)  U(4) 

READ(16,')  U(5) 

READ(15,')  U(6) 

READ(15,*)  U(7) 

READ(15,*)  VTRFPS 
U(8)=VTRFPS/1000.0cl0 


C 

C  PAR(1)=DELESD 

C  PAR(2)=DRUDD  THE  PARAMETERS,  IN  DEGREES 
C  PAR(3)=DDA 

C  PAR(4)=ENGPA  PORT  ENGINE  THRUST.  POUNDS/1 000 
C  PAR(5)=ENGSA  STARBORD  ENGINE  THRUST,  POUNDS/1 000 
C  PAR(6)=TPTAL  PITCH  THRUST  VECTOR,  DEG 
C  PAR(7)=TYTAL  YAW  THRUST  VECTOR,  DEG 
C  PAR(8)=TTHRST  TOTAL  THRUST,  POUNDS/1 000 

C  MODIFIED  13  AUG  91 
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G  PAR(9)=DELEDD  DIFFERENTIAL  ELEVATOR,  IN  DEGREES 
C  PAR(10)=BLEND  TRANSITION  PARAMETER  FROM  STATIC  TO 
C  ROTARY  BAUNCE  COEFICIENT  DATA 

C 

READ(15,*)  PAR(1) 

READ(15,*)  PAR(2) 

READ{15.*)  PAR(3) 

READ{15,')  PAR(4) 

READ(I5,*)  PAR(5) 

READ(15.*)  PAR(6) 

READ(15.')  PAR(7) 

READ{15,*)  PAR(8) 

C  MODIFIED  13  AUG  91 

READ(15.*)  PAR(9) 

READ(15.')  PAR(IO) 


RETURN 

END 


SUBROUTINE  INIT 


IMPLICIT  DOUBLE  PRECISIGN(A-H.O-2) 

C 

COMMON  /BLCSS/NDIM,ITMX,NPAR,ICP,IID,NMX,IPS,IRS 
COMMON  /BLCPS/  NTST.NCOL.IANCH,NMXPS,IAD,NPR,NWTN,ISP,ISWl 
COMMON  /BLDLS/  DS,DSMIN,DSMAX,IADS 
COMMON  /BLUM/  RL0,RL1,AQ,A1,PAR{10) 

COMMON  /BLOPT/  ITNW,MXBF.IPLT,ICP2,ILP 
COMMON  /BLEPS/  EPSU.EPSL.EPSS.EPSR 
C 
C 

C  IN  THIS  SUBROUTINE  THE  USER  SHOULD  SET  THOSE  CONSTANTS 
C  THAT  REQUIRE  VALUES  DIFFERENT  FROM  THE  DEFAULT  VALUES 
C  ASSIGNED  IN  THE  LIBRARY  SUBROUTINE  DFINIT.  FOR  A 
C  DESCRIPTION  OF  THESE  CONSTANTS  SEE  THE  DOCUMENTATION 
C  CONTAINED  IN  THE  LIBRARY.  COMMON  BLOCKS  CORRESPONDING  TO 
C  CONSTANTS  THAT  THE  USER  WANTS  TO  CHANGE  MUST  BE  INSERTED 
C  ABOVE.  THESE  COMMON  BLOCKS  SHOULD  OF  COURSE  BE  IDENTICAL 
C  TO  THOSE  IN  DFINIT. 

C 

C 

DSMAX  =  lO.OdO 
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DSMIN  =0.0000001  OdO 
EPSU  =  1  .OD-07 
EPSL  =  1.0D-07 
EPSS  =  1.0D-05 
EPSR  =1.0D-07 
lAD  =  1 
ILP  =  1 
ITMX  =40 
ITNW  =20 
MXBF  =5 
NDIM  =8 
NPAR  =  10 

OPEN(UNIT=25,FILE=’fort.25’) 

REWIND  (25) 

READ(25.*)  RL0.RL1 
READ{25,*)  A0.A1 
READ{25.*)  DS 
READ{25,*)  NMX 

READ{25,*)  NTST.NCOL.NMXPS.NPR 

READ(25.‘)  ISP.IRS.ICP.ICP2.IPLT,IPS 

READ(25,*)  ISW1 

RETURN 

END 

SUBROUTINE  BOND 


RETURN 

END 

SUBROUTINE  ICND 


RETURN 

END 


SUBROUTINE  COEFF(U,PAR,NDIM.ICP) 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  /ACDATA/  BWING.CWING.SREF.RHO.RMASS 
COMMON  /SEIZE/  CX,CY,CZ,CLM,CMM.CNM 
COMMON /SEIZET/CXT,CYT,CZT,CLMT,CMMT,CNMT 
DIMENSION  U(NDIM),PAR(10) 


THE  PRIMARY  SOURCE  OF  THESE  COEFFICIENT  EQUATIONS  IS 
SUBROUTINE  AROlO  FROM  MCAIR  CODE  USED  IN  THE  F15 
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C  BASELINE  SIMULATOR. 

C 

C  MOST  OF  THE  COEFFICIENTS  USED  IN  THE  EQUATIONS  WERE 
C  COMPUTED  USING  SAS  WITH  RAW  DATA  FROM  THE  FI 5  SIMULATOR 
C  DATA  TABLES. 

C 

C  THIS  SUBROUTINE  IS  CALLED  BY  THE  DRIVER  PROGRAM  FOR  THE 
C  AUTO  SOFTWARE.  IT  MERELY  TAKES  INPUTS  ON  THE  A/C 
C  STATE,  CONTROL  SURFACE  POSITIONS.  AND  THRUST  VALUES 
C  AND  RETURNS  THE  APPROPRIATE  AERO  COEFFICIENTS  CX,  CY, 

C  CZ,  CL,  CM,  AND  CN. 

C 

C  INPUTS  TO  THIS  SUBROUTINE 
C 

C  AL  -  ANGLE  OF  ATTACK,  DEG 
C  BETA  -  SIDESLIP  ANGLE,  DEG 
C  DDA  -  AILERON  DEFLECTION  ANGLE,  DEG 
C  DELEDD  -  DIFFERENTIAL  TAIL  DEFLECTION  ANGLE,  DEG 
C  DELESD  -  SYMMETRICAL  TAIL  DEFLECTION  ANGLE,  DEG 
C  DRUDD  -  RUDDER  DEFLECTION,  POSITIVE  TRAILING  EDGE 
C  LEFT.  DEG 

C  P  -  ROLL  RATE.  RAD/SEC 

C  Q  -  PITCH  RATE,  RAD/SEC 

C  R  -  YAW  RATE,  RAD/SEC 

C  ENGPA  -PORT  ENGINE  THRUST,  POUNDS/1 000 
C  ENGSA  -STARBOARD  ENGINE  THRUST,  POUNDS/1 000 
C  TYTAL  -  YAW  THRUST  ANGLE,  DEG 
C  TPTAL  -  PITCH  THRUST  ANGLE.  DEG 
C  TTHRST  -  TOTAL  THRUST,  POUNDS/1000 
C  VTRFPS  -  TRUE  AIRSPEED,  FT/SEC 
C 

C  INTERMEDIATE  VARIABLES  USED  IN  THIS  SUBROUTINE 
C 

C  ABET  -  ABSOLUTE  VALUE  OF  BETA,  DEG 
C  ARUD  -  ABSOLUTE  VALUE  OF  RUDDER  DEFLECTION,  DEG 
C  SWING  -  WING  SPAN,  FEET 
C  CA  -  COSINE  RAL  (RAL  IN  RADIANS) 

C  CD  -  COEFFICIENT  OF  DRAG 
C  CL  -  BASIC  LIFT  COEFFICIENT 
C  CWING  -  MEAN  AERODYNAMIC  CHORD,  FEET 
C  DAHD  -  DIFFERENTIAL  ELEVATOR  DEFLECTION.  DEG 
C  DAHLD  -  LEFT  AILERON  DEFLECTION,  DEG 
C  DAHRD  -  RIGHT  AILERON  DELFECTION,  DEG 
C  DELEDR  -  DIFFERENTIAL  TAIL  DEFLECTION  ANGLE,  RAD 
C  DELESR  -  SYMMETRIC  TAIL  DEFLECTION  ANGLE,  RAD 
C  ENGP  -  PORT  ENGINE  THRUST,  POUNDS 
C  ENGS  -  STARBOARD  ENGINE  THRUST,  POUNDS 
C  PTAL  -  PITCH  THRUST  VECTOR,  RAD 

C  QBARS  -  DYNAMIC  PRESSURE  TIMES  WING  REFERENCE  AREA. 

C  LBF 
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RABET  -  ABSOLUTE  VALUE  OF  BETA,  RADIANS 
RAL  -  ABSOLUTE  VALUE  OF  ALPHA.  RADIANS 
RARUD  -  ABSOLUTE  VALUE  OF  RUDDER,  RADIANS 
SA  -  SINE  RAL  (RAL  IN  RADIANS) 

YTAL  -  YAW  THRUST  VECTOR,  RAD 


C  OUTPUTS  FROM  THIS  SUBROUTINE 
C 

C  CX  -  BASIC  AXIAL  FORCE  COEFFICIENT,  BODY  AXIS, 

C  +  FORWARD 

C  CY  -  BASIC  SIDE  FORCE  COEFFICIENT,  BODY  AXIS, 

C  +  RIGHT 

C  CZ  -  BASIC  NORMAL  FORCE  COEFFICIENT,  BODY  AXIS, 

C  +  DOWN 

C  CLM  -  BASIC  ROLLING  MOMENT  COEFFICIENT,  BODY  AXIS, 

C  +  R  WING  DOWN 

C  CMM  -  BASIC  PITCHING  MOMENT  COEFFICIENT,  BODY  AXIS, 

C  +  NOSE  UP 

C  CNM  -  BASIC  YAWING  MOMENT  COEFFICIENT,  BODY  AXIS, 

C  +  NOSE  RIGHT 

C 

C  ANGLES  USED  IN  CALCULATING  CL,  CLLDB, ...,  ARE  IN 
C  RADIANS.  THIS  IS  BECAUSE  RADIANS  WERE  USED  IN  THE  CURVE 
C  FITTING  PROGRAM  TO  OBTAIN  THE  COEFFICIENTS  OF  THE 

C  ALPHA.  BETA . TERMS  IN  THE 

C  FOLLOWING  EQUATIONS. 

C 

C 

C  MOMENT  REFERENCE  CENTER  WAS  SET  IN  AROlO  PROGRAM  AS: 

C 

C  DATA  CMCGR  /.2565/.  CNCGR  /.2565/ 

C 

C  THE  AERO  STABILITY  DATA  WAS  TAKEN  REFERENCED  TO  THESE  CG 
C  LOCATIONS.  THE  MOMENTS  OF  INERTIA  AND  OTHER  AIRCRAFT 
C  DATA  ARE  FOR  A  CLEAN  CONFIGURATION  TEST  AIRCRAFT  WITH  A 
C  CG  AT  THE  SAME  CG.  AS  A  RESULT.  THERE  IS  NO  ’CG  OFFSET’ 

C  TO  BE  COMPUTED. 

C 

IWRITE=0 

C 

AL=U(1) 

BETA=U{2) 

P=U(3) 

Q=U(4) 

R=U(5} 

THETA=U(6) 

PHI=U(7) 
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VTRFPS=U(8)*1000. 

c 

DELESD=PAR(1) 

DRUDD=PAR(2) 

DDA=PAR(3) 

ENGPA=PAR(4) 

ENGSA=PAR(5) 

TPTAL=PAR(6) 

TYTAL=PAR(7) 

TTHRST=PAR(8) 

C 

DEGRAD=57.29577951 
DELESR=DELESD/DEGRAD 
YTAL=TYTAL/DEGRAD 
PTAL=TPTAUDEGRAD 

IF  BLOCK  TO  CHANGE  TOTAL  THRUST 

IF(ICP.EQ.8)THEN 
DIFT=PAR(4)-PAR(5) 

THALF=TTHRST/2.0d0 
ENGPA=THALf+DIFT/2.0d0 
ENGSA=THALF-DIFT/2.0d0 
ENDIF 

ENGP=ENGPA*1 000.0 
ENGS=ENGSA'1 000.0 

QBARS=0.5dO*RHO*VTRFPS*VTRFPS‘SREF 
CO2V=CWING/(2.0d0*VTRFPS) 
BO2V=BWING/(2.0d0*VTRFPS) 
QSB=BWING*QBARS 
ARUD=ABS(DRUDD) 
RARUD=ARUD/DEGRAD 
RAL=AL/DEGRAD 
ABET=ABS(BETA) 

RABET=ABET/DEGRAD 


1)  ALL  THE  AERODYNAMIC  COEFFICIENTS  IN 
THIS  VERSION  OF  THE  DRIVER  PROGRAM 
ARE  TAKEN  DIRECTLY  FROM  THE  1988 
FI  5  AEROBASE  (0.6  MACH,  20000  FEET) 
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2)  THIS  SECTION  SUMMARIZES  THE 
AERODYNAMIC  COEFFICIENTS  AS  TO  WHAT 
THEY  ARE  AND  HOW  THEY  ARE  USED. 

THE  FIRST  ACCRONYM  IS  THE  JOVIAL 
NAME  OF  THE  AERODYNAMIC  COEFFICIENT 
(CFX1.  ETC),  THE  SECOND  ACCRONYM  IS 
THE  F15  AEROBASE  CODE  OR  CTAB  NAME 
(ATAB15.  ETC).  A  BRIEF  DEFINITION 

OF  THE  AERODYNAMIC  COEFFICIENT  IS 
ALSO  PROVIDED. 

3)  THERE  IS  ALSO  A  SECTION  THAT 
PROVIDES  A  TABLE  OF  CONVERSIONS 
BETWEEN  WHAT  THE  VARIABLE  IS  CALLED 
IN  THE  OR  SECTION  OF  THIS  PROGRAM 
AND  ITS  NAME  IN  THE  1988  FI 5 
AEROBASE.  FOR.THE  SAKE  OF 
CONTINUITY  THE  ORIGINAL  PROGRAM  NAME 
IS  USED  AND  THE  1988  F15  AEROBASE 
NAME  IS  PROVIDED  AS  BOOK  KEEPING 
INFORMATION. 


***•#»*****♦»***#*♦**♦♦***#•**********» 


CFX  =  FORCE  IN  STABILITY  AXIS  X  DIRECTION  (CD  IN  BODY  AXIS) 
(FUNCTION  OF  CL  OR  CFZI) 

CFX  =  CFX1  +  CXRB  +  STORE  INCREMENTS  +  CXDSPD  +  DCXLG  +  DCD 

CFXI  =  ATAB15  =  PERFORMANCE  DRAG  COEFFICIENT  -  CD 
CXRB  =  ATAB22  =  DELTA  CD  DUE  TO  CG  (=0.0) 

CXDSPD  =  ATAB27  =  DELTA  CD  DUE  TO  SPEEDBRAKE  (NORMALLY  =  0.0436) 
SET  TO  0  SINCE  THIS  STUDY  IS  CONCERNED 
WITH  HIGH  ANGLES 

OF  ATTACK  PHENOMENON  (>40  DEGREES)  AND  BECAUSE 
THE  SPEEDBRAKE  WILL  NOT  DEPLOY  AT  ANGLES  OF 
ATTACK  GREATER  THAN  15  DEGREES. 

DCXLG  =  ATAB19  =  DELTA  CD  DUE  TO  REYNOLD’S  NUMBER  (=-0.0005) 

DCD  =  BTAB03  =  DELTA  CD  DUE  TO  2-PLACE  CANOPY  (F15B)  (=0.0005) 

* . NOTE  THAT  DCXLG  AND  DCD  CANCEL  EACH  OTHER  *“**** 


CFY  =  FORCE  IN  BODY  AXIS  Y  DIRECTION 

CFY  =  CFY1*EPA02  +  CYDAD'DAILD  +  (CYDRD*DRUDD*DRFLX5]*EPA43 
+(CYDTD*DTFLX5  +  DTFLX6)*DTALD  +  CFYP*PB  +  CFYR‘RB 
+CYRB  +  STORE  INCREMENTS  +  DCYB'BETA 
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CFY1  =  ATAB16  =  BASIC  SIDE  FORCE  COEFFICIENT  -  CY(BETA) 

EPA02  =  ATAB21  =  BETA  MULTIPLIER  TABLE 
CYDAD  =  ATAB75  =  SIDE  FORCE  COEFFICIENT  DUE  TO  AILERON 
DEFLECTION 

DAILD  =  AILERON  DEFLECTION  (DEG) 

CYDRD  =  ATAB69  =  SIDE  FORCE  COEFFICIENT  DUE  TO  RUDDER  DEFLECTION 
DRUDD  =  RUDDER  DEFLECTION  (DEG) 

DRFLX5  =  ATAB88  =  FLEX  MULTIPLIER  ON  CYDRD  (=0.89) 

EPA43  =  ATAB30  =  MULTIPLIER  ON  CNDR,  CLDR,  CYDR  DUE  TO 
SPEEDBRAKE 
(=1.0) 

CYDTD  =  ATAB72  =  SIDE  FORCE  COEFFICIENT  DUE  TO  DIFFERENTIAL  TAIL 
DEFLECTION  -  CYDDT 

DTFLX5  =  ATAB10  =  FLEX  MULTIPLIER  ON  CYDTD  (=0.975) 

DTFLX6  =  ATAB77  =  FLEX  INCREMENT  TO  CYDTD  (=0.0) 

DTALD  =  DIFFERENTIAL  TAIL  DEFLECTION  (DEG)  WHICH  IS 
DIRECTLY  PROPORTIONAL  TO  AILERON  DEFLECTION 
AND  IS  PRIMARILY  USED  TO  ASSIST  IN  ROLLING  THE 
F-15B  (DTALD=0.3*DAILD) 

CFYP  =  ATAB13  =  SIDE  FORCE  COEFFICIENT  DUE  TO  ROLL  RATE  (CYP) 

PB  =  (PEOBB*SPAN)/(2*VILWF) 

PEOBB  =  ROLL  RATE  IN  RAD/SEC  =  P 
SPAN  =  WING  SPAN  =  42.8  FEET  =  BWING 
VILWF  =  VELOCITY  IN  FT/SEC  =  VTRFPS 
CFYR  =  ATAB07  =  SIDE  FORCE  COEFFICIENT  DUE  TO  YAW  RATE  (CYR) 

RB  =  (REOBB*SPAN)/(2*VILWF) 

REOBB  =  YAW  RATE  IN  RAD/SEC  =  R 
CYRB  =  ATAB93  =  ASSYMETRIC  CY  AT  HIGH  ALPHA  (ANGLE  OF  ATTACK) 
DCYB  =  0.0  THERE  IS  NO  INCREMENT  DELTA  CYB  (SIDE 
FORCE) 

DUE  TO  A  2-PLACE  CANOPY  ON  THE  F15B.  THIS  IS 
BECAUSE  THE  SAME  CANOPY  IS  USED  ON  BOTH  THE 
BASELINE  F15A  AND  THE  F15B.  THE  SIDEFORCE  IS 
THE  SAME  FOR  BOTH  VERSIONS  OF  THE  F15  AND 
ALREADY  INCLUDED  IN  THE  BASIC  SIDE  FORCE  (CFYI). 

THE  TWO  PLACE  CANOPY  IS  MOUNTED  DIFFERENTLY 
HOWEVER,  SO  THERE  IS  A  DIFFERENCE  IN  YAWING  AND 
C  ROLLING  MOMENT. 

C  (SEE  DCNB  AND  DCLB) 

C 

C 

C 

c 

C  CFZ  =  FORCE  IN  STABILITY  AXIS  Z  DIRECTION  (CL  IN  BODY  AXIS) 

C  CFZ  =  CFZ1  +  CZDSPD  +  STORE  INCREMENTS  +  DCL'BETA 

C 

C 

C  CFZ1  =  ATAB17  =  BASIC  LIFT  COEFFICIENT  -  CL 
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CZDSPD  =  ATAB26  =  DELTA  CL  DUE  TO  SPEEDBRAKE 

SET  TO  0  DUE  TO  THE  REASONS  GIVEN  ABOVE  IN 
CXDSPD 

DCL  =  BTAB01  =  DELTA  CL  DUE  TO  2-PLACE  CANOPY  (F15B)  (=0.0) 


It  It**************  **********  It  tit*  1i*********ltlt***it1e***i,*-k**lr1e1tltltint1Ht* 


CML  =  TOTAL  ROLLING  MOMENT  COEFFICIENT  IN  BODY  AXIS 
CML  =  CML1‘EPA02  +  CLDAD*DAILD  +  (CLDRD*DRUDD*DRFLX1)*EPA43  + 
(CLDTD'DTFLXI  +  DTFLX2)‘DTALD  +  CMLP*PB  +  CMLR'RB  + 

STORE  INCREMENTS  +  CLDSPD  +  DCLB*BETA 


C 

C  CML1  =  ATAB01  =  BASIC  ROLLING  MOMENT  COEFFICIENT  -  CL(BETA) 

C  EPA02  =  ATAB21  =  BETA  MULTIPLIER  TABLE 
C  CLDAD  =  ATAB73  =  ROLL  MOMENT  COEFFICIENT  DUE  TO  AILERON 
C  DEFLECTION 

C  -(CLDA) 

C  DAILD  =  AILERON  DEFLECTION  (DEG) 

C  CLDRD  =  ATAB67  =  ROLLING  MOMENT  COEFFICIENT  DUE  TO  RUDDER 
C  DEFLECTION  -(CLD) 

C  DRUDD  =  RUDDER  DEFLECTION  (DEG) 

C  DRFLXI  =  ATAB80  =  FLEX  MULTIPLIER  ON  CLDRD  (=0.85) 

C  EPA43  =  ATAB30  =  MULTIPLIER  ON  CNDR.  CLDR,  CYDR  DUE  TO 
C  SPEEDBRAKE 

C  (=1.0) 

C  CLDTD  =  ATAB70  =  ROLL  MOMENT  COEFFICIENT  DUE  TO  DIFFERENTIAL 
C  TAIL 

C  DEFLECTION  -  CLDD 

C  DTFLX1  =  ATAB04  =  FLEX  MULTIPLIER  ON  CLDTD  (=0.975) 

C  DTFLX2  =  ATAB84  =  FLEX  INCREMENT  TO  CLDTD  (=0.0) 

C  DTALD  =  DIFFERENTIAL  TAIL  DEFLECTION  (DEG)  WHICH  IS 
C  DIRECTLY  PROPORTIONAL  TO  AILERON  DEFLECTION 

C  AND  IS  PRIMARILY  USED  TO  ASSIST  IN  ROLLING  THE 

C  F-15B 

C  (DTALD  =  0.3*DAILD) 

C  CMLP  =  ATAB02  =  ROLL  DAMPING  DERIVATIVE  -CLP 
C  PB  =  (PEOBB*SPAN)/(2'VILWF) 

C  PEOBB  =  ROLL  RATE  IN  RAD/SEC  =  P 

C  SPAN  =  WING  SPAN  =  42.8  FEET  =  BWING 

C  VILWF  =  VELOCITY  IN  FT/SEC  =  VTRFPS 

C  CMLR  =  ATAB11  =  ROLLING  MOMENT  COEFFICIENT  DUE  TO  YAW  RATE  -  CLR 
C  RB  =  (REOBB*SPAN)/(2'VILWF) 

C  REOBB  =  YAW  RATE  IN  RAD/SEC  =  R 

C  CLDSPD  =  ATAB29  =  DELTA  CL  DUE  TO  SPEEDBRAKE 
C  SET  TO  0  DUE  TO  THE  REASONS  GIVEN  ABOVE  IN 

C  CXDSPD 

C  DCLB  =  BTAB04  =  INCREMENT  DELTA  CLB  (ROLLING  MOMENT)  DUE  TO 
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2-PLACE 

CANOPY  FROM  PSWT  499 


****************************************************************** 


CMM  =  TOTAL  PITCHING  MOMENT  COEFFICIENT  IN  STABILITY  AXIS 
(BODY  AXIS- AS -WELL) 

CMM  =  CMM1  +  CMMQ'QB  +  STORE  INCREMENTS  +  CMDSPD  +  DCM 

CMM1  =  ATAB03  =  BASIC  PITCHING  MOMENT  COEFFICIENT  -  CM 
CMMQ  =  ATAB05  =  PITCH  DAMPING  DERIVATIVE  -  CMQ 
OB  =  (QEOBB*MAC)/(2'VILWF) 

QEOBB  =  PITCH  RATE  IN  RAD/SEC  =  Q 
C  MAC  =  MEAN  AERODYNAMIC  CHORD  =  15.94  FEET  =  OWING 

C  VILWF  =  VELOCITY  IN  FT/SEC  =  VTRFPS 

C  CMDSPD  =  ATAB25  =  DELTA  CM  DUE  TO  SPEEDBRAKE 
C  SET  TO  0  DUE  THE  REASONS  GIVEN  ABOVE  IN  CXDSPD 

C  DCM  BTAB02  =  DELTA  CM  DUE  TO  2-PLACE  CANOPY  (F15B)  (=0.0) 

C 

C 

Q** ********** *********************************************** ********** 

C 

c 

c 

C  CMN  =  TOTAL  YAWING  MOMENT  COEFFICIENT  IN  BODY  AXIS 
C  CMN  =  CMNrEPA02  +  CNDAD'DAILD  +  (CNDRD*DRUDD‘DRFLX3]*EPA43 
C  +(CNDTD*DTLX3  +  DTFLX4]*DTALD  +  CMNP*PB  +  CMNR'RB  +  CNRB 
C  +DCNB2*EPA36  +  STORE  INCREMENTS  +  CNDSPD  +  DCNB‘BETA 
C 
C 

C  CMN1  =  ATAB12  =  BASIC  YAWING  MOMENT  COEFFICIENT  -  CN  (BETA) 

C  EPA02  =  ATAB21  =  BETA  MULTIPLIER  TABLE 
C  CNDAD  =  ATAB74  =  YAW  MOMENT  COEFFICIENT  DUE  TO  AILERON 
C  DEFLECTION  -CNDA 

C  DAILD=  =  AILERON  DEFLECTION  (DEG) 

C  CNDRD  =  ATAB68  =  YAWING  MOMENT  COEFFICIENT  DUE  TO  RUDDER 
C  DEFLECTION  -CNDR 

C  DRUDD  =  RUDDER  DEFLECTION  (DEG) 

C  DRFLX3  =  ATAB85  =  FLEX  MULTIPLIER  ON  CNDRD 

C  EPA43  =  ATAB30  =  MULTIPLIER  ON  CNDR,  CLDR,  CYDR  DUE  TO  SPEEDBRAKE 
C  CNDTD  =  ATAB71  =  YAWING  MOMENT  COEFFICIENT  DUE  TO  DIFFERENTIAL  TAIL 
C  DEFLECTION  -  CNDDT 

C  DTFLX3  =  ATAB08  =  FLEX  MULTIPLIER  ON  CNDTD 
C  DTFLX4  =  ATAB09  =  FLEX  INCREMENT  ON  CNDTD  (=0.0) 

C  DTALD  =  =  DIFFERENTIAL  TAIL  DEFLECTION  (DEG)  WHICH  IS 

C  DIRECTLY  PROPORTIONAL  TO  AILERON  DEFLECTION 

C  AND  IS  PRIMARILY  USED  TO  ASSIST  IN  ROLLING 

C  THE  F-1 5B  (DTALD  =  0.3*DAILD) 
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C  CMNP  =  ATAB06  =  YAWING  MOMENT  COEFFICIENT  DUE  TO  ROLL  RATE  -  CNF 
C  PB  =  (PEOBB*SPAN)/(2*VILWF) 

C  PEOBB=ROLL  RATE  IN  RAD/SEC  =  P 

C  SPAN  =  WING  SPAN  =  42.8  FT  =  BWING 

C  VILWF  =  VELOCITY  IN  FT/SEC  =  VTRFPS 

C  CMNR  =  ATAB14  =  YAW- DAMPING  DERIVATIVE  -  CNR 
C  RB  =  (REOBB*SPAN)/(2*VILWF) 

C  REOBB  =  YAW  RATE  IN  RAD/SEC  =  R 

C  CNRB  =  ATAB86  =  ASSYMETRIC  CN  AT  HIGH  ALPHA 

C  DCNB2  =  ATAB44  =  DELTA  CNB  WITH  STABILATOR  EFFECT  -  DELCNB  (=0.0) 

C  EPA36  =  ATAB94  =  MULTIPLIER  ON  DCNB2  (=BETA) 

C  CNDSPD  =  ATAB28  =  DELTA  CN  DUE  TO  SPEEDBRAKE 
C  SET  TO  0  DUE  TO  THE  REASONS  GIVEN  ABOVE  IN  CXDSPD 

C  DCNB  =  BTAB05  =  INCREMENT  DELTA  CNB  (YAWING  MOMENT)  DUE  TO  . 

C  2-PLACE  CANOPY  (FI  5B) 

C 

C 

^»************»****»#******i»i»**#*******i»****»*****»»»********  ********* 

C 

c 

c 

C  MISCELLANEOUS  COEFFICIENTS  AND  NAME  CONVERSION  TABLE 


C 

.C 


c 

1988  FI  5 

ORIGINAL 

c 

c 

AEROBASE  NAME 
************* 

PROGRAM  NAME  DEFINITION 

************  •»♦**♦**** 

c 

c 

AL77D 

AL 

ANGLE  OF  ATTACK 

c 

(DEG) 

c 

BE77D 

BETA 

SIDESLIP  ANGLE 

c 

(DEG) 

c 

BE77D 

RBETA 

SIDESLIP  ANGLE 

c 

(RAD) 

c 

B077D 

ABET 

ABSOLUTE  VALUE  OF 

c 

SIDESLIP  ANGLE 

c 

(DEG) 

c 

DAILA 

DAILA 

ABSOLUTE  VALUE  OF 

c 

AILERON  DEFLEC¬ 

c 

TION  (DEG) 

c 

DAILD 

DDA 

AILERON  DEFLEC¬ 

c 

TION  (DEG) 

c 

DRUABS 

ARUD 

ABSOLUTE  VALUE  OF 

c 

RUDDER  DEFLEC¬ 

c 

TION  (DEG) 

c 

DRUABS 

RARUD 

ABSOLUTE  VALUE  OF 

c 

RUDDER  DEFLEC¬ 

c 

TION  (RAb) 

c 

DRUDD 

DRUDD 

RUDDER  DEFLECTION 

c 

(DEG) 
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DSTBD  DELESD(R)  AVERAGE 

STABILATOR 
DEFLECTION 
DEG  (RAD) 

DTALD  DELEDD(R)  DIFFERENTIAL  TAIL 

DEFLECTION 
DEG  (RAD) 


RBETA=BETA/DEGRAD 

DAILA=ABS(DDA) 


PB=(P'BWING)/(2.0d0*VTRFPS) 

QB=(Q*CWING)/(2.0ci0*VTRFPS) 

RB=(R*BWING)/(2.0ciO*VTRFPS) 

C 

C  THE  F-15B  AERO  DATA  TABLES  DO  NOT  CONTAIN  STABILITY  COEFFICIENT 
C  DATA  FOR  BETA  AND  RUDDER  DEFLECTION  ,DRUDD,  LESS  THAN  0 
C  DEGREES.  THE  ABSOLUTE  VALUE  OF  BETA,  ABET,  AND  THE  ABSOLUTE 
C  VALUE  OF  RUDDER  DEFLECTION,  ARUDD,  ARE  USED  IN  THE  FOLLOWING 
C  EG'  ATIONS.  IN  RADIANS  THESE  PARAMETERS  ARE  RABET  AND  RARUD, 

C  RESPECTIVELY.  IN  SOME  CASES  THE  COEFFICIENT  IS  MULTIPLIED  BY  A 
C  -1  FOR  PARAMETER  VALUES  LESS  THAN  ZERO. 

C 

C  EPA02  IS  A  MULTIPLIER  THAT  ADJUSTS  THE  PARTICULAR  COEFFICIENT 
C  IT  IS  WORKING  ON  (CFY1  .CML1  ,CMN1)  BY  CHANGING  THAT  PARTICULAR 
C  COEFFICIENTS  SIGN  (POSITIVE  OR  NEGATIVE)  DEPENDENT  ON  THE  SIGN 
C  OF  THE  SIDESLIP  ANGLE  (BETA).  IF  BETA  IS  NEGATIVE  THEN 
C  EPA02=-1.0.  IF  BETA  IS  POSITIVE  THEN  EPA02=1.0.  SINCE  THIS 
C  FUNCTION  IS  DISCONTINUOUS  AT  THE  ORIGIN  A  CUBIC  SPLINE  HAS 
C  BEEN  EMPLOYED  TO  REPRESENT  THIS  FUNCTION  IN  ORDER  THAT 
C  AUTO  CAN  RUN. 

C 

C 

IF  (BETA  .LT.  -1.0)  THEN 
EPA02S=  -I.OdO 
ENDIF 
C 

IF  ((BETA  .GE.  -1.0)  .AND.  (BETA  .LE.  1.0))  THEN 
EPA02S=-1 .0d0+(1 .50d0*((BETA+1 .0d0)'*2.0d0))- 
1  (0.50d0*((BETA+1.0d0)*'3.0d0)) 

ENDIF 

C 

IF  (BETA  .GT.  1.0)  THEN 
EPA02S=1  .OdO 
ENDIF 
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IF  (BETA  .LT.  -5.0)  THEN 
EPA02L= -I.OdO 
ENDIF 
C 

IF  ((BETA  .GE.  -5.0)  .AND.  (BETA  .LE.  5.0))  THEN 
EPA02L=-1.0cl0+(0.060cl0‘((BETA+5.0d0)**2.0d0))- 
1  (O.O04OdO*((BETA+5.OdO)“3.OdO)) 

ENDIF 

C 

IF  (BETA  .GT.  5.0)  THEN 

EPA02L=1.0d0 

ENDIF 


DTALD=0.30d0’DAILD 

DELEDD=0.30d0*DDA 

DELEDR=0.30d0*(DDA/DEGRAD) 


CFZ1=-0.00369376+(3.78028702*RAL)+(0.6921459*RAL*RAL)-(5.0005867 
+*(RAL**3))+(1 .944781 99*(RAL"4))+(0.40781 955'DELESR)+(0.1 01 1 4579 
+*(DELESR*DELESR)) 

CF2=CFZ1 


CL=CFZ1 757.29578 

THIS  CONVERSION  OF  CFZ1  TO  CL  IS  AN  ARTIFACT  FROM  THE 
CURVE  FITTING  PROCESS  WHERE  ALL  THE  INDEPENDENT  VARIABLES 
WERE  ANGLES  THAT  WERE  CONVERTED  FROM  DEGREES  TO  RADIANS. 
IT  JUST  SO  HAPPENED  THAT  FOR  CFXl  ONE  OF  THE  VARIABLES 
WAS  NOT  AN  ANGLE  BUT  A  DIMENSIONLESS  COEFFICIENT. 


CFXl  =0.01 806821 +(0.01 556573*CL)+(498.96208868*CL*CL) 
+-(1 4451 .5651 8396*(CL**3))+(21 32344.61 84755*(CL"4)) 
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TRANSITIONING  FROf\^  LOW  AOA  DRAG  TABLE  TO  HIGH  AOA  DRAG  TABLE 

CFX2=0.0267297-{0.1 064691 9*RAL)+(5.39836337*RAL*RAL) 
+-(5.0086893*(RAL**3))+(1 .341 481 93*(RAL**4))+ 
+(0.20978902*DELESR)+(0.30604211*(DELESR**2))+0.09833617 

A1=20.0d0/DEGRAD 

A2=30.0ci0/DEGRAD 

A12=A1+A2 

BA=2.0/(-Ar*3+3.*Al*A2*(Al-A2)+A2**3) 

BB=-3.0d0*BA*(A1  +A2)/2.0d0 
BC=3.0d0*BA*A1*A2 
BD=BA*A2**2*(A2-3.0d0*A1  )/2.0d0 
F1=BA*RAL**3+BB*RAL*'2+EC’'RAL+BD 
F2=-BA*RAL**3+(3.0dO*Al2*BA+BB)’RAL**2.0dO- 
+  (BC+2.0d0*A1 2*BB+3.0d0*A1 2**2*BA)*RAL+ 

+  BD+A12*BC+A12'*2*BB+A12**3*BA 

IF  (RAL  .LT.  A1)THEN 

CFX=CFX1 

ELSEIF  (RAL  .GT.  A2)  THEN 
CFX=CFX2 
ELSE 

CFX=CFXrF1.+CFX2*F2 

ENDIF 


DTFLX5=0.975d0 

DRFLX5=0.89d0 

C 

CFY1  =-0.05060386-(0.1 2342073*RAL)+{1 .04501 1 36*RAL*  RAL) 

+-(0.1 723951 6*(RAL*’3))-(2.90979277'(RAL**4)) 
++(3.06782935*(RAL"5))-(0.884221 1 6*(RAL**6)) 
+-(0.06578812'RAL*RABET)-(0.71521988*RABET)-(0.00000475273 
+*(RABET**2))-(0.04856168*RAL*DELESR)-(0.05943607*RABErDELESR)+ 
+(0.0201 8534*DELESR) 

C 

IF  (RAL  .LT.  .52359998)  THEN 
C 

CFYP=0.014606188+(2.52405055'RAL)-(5.02687473*(RAL**2)) 

+-(106.43222962'(RAL"3))+(256.80215423*(RAL*'4)) 
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++(1 256.39636248‘(RAL**5)) 

+-(3887.928781 73*(RAL**6))-(2863.1 6083460‘{RAL**7))+ 

+(1 7382.72226362*(RAL**8))-(1 3731 .65408408*(RAL**9)) 

ENDIF 

C 

IF  ((RAL  .GE.  .52359998)  .AND.  (RAL  .LE.  .610865))  THEN 
C 

CFYP=0.0023651 1  +(0.52044678*(RAL-0.52359998))-(1 2.8597002*(RAL- 
+0.52359998)**2)+(75.46138*(RAL-0.52359998)**3) 

ENDIF 

C 

IF  (RAL  .GT.  0.610865)  THEN 
C 

CFYP=0.0d0 

ENDIF 

C 

IF  (RAL  .LT.  -0.06981)  THEN 
C 

CFYR=0.35d0 

ENDIF 

C 

IF  ((RAL  .GE.  -0.06981)  .AND.  (RAL  .LT,  0.0))  THEN 
C 

CFYR=0.34999999+(35.4012413*(RAL+0.06981)**2)-(493.33441162* 

+(RAL+0.06981)**3) 

ENDIF 

C 

IF  ((RAL  .GE.  0.0)  .AND.  (RAL  .LE.  0.523599))  THEN 
C 

CFYR=0.35468605-(2.26998141*RAL)+(51.82178387'RAL'RAL) 

+-(71 8.55069823*(RAL"3)) 

++(4570.004921 72*(RAL*'4))-(1 4471 .8802835r(RAL*’5))+ 
+(22026.58930662'(RAL**6))-(127S5.99029404*{RAL**7)) 

ENDIF 

C 

IF  ((RAL  .GT.  0.523599)  .AND.  (RAL  .LE.  0.61087))  THEN 
C 

CFYR=0.00193787+(1.78332496*(RAL-0.52359903))-(41.63l98853*(RAL- 

+0.52359903)**2)+(239.97909546'(RAL-0.52359903)"3) 

ENDIF 

C 

IF  (RAL  .GT.  0.61087)  THEN 
C 

CFYR=0.0d0 

ENDIF 

C 

IF  (RAL  .LT.  0.55851)  THEN 
C 

CYDAD=-0.00020812+(0.00062122*RAL)+{0.00260729*RAL'RAL) 
++(0.00745739*(RAL**3))-(0.036561 1  *(RAL**4)) 


131 


oooooooooooo 


+-(0.04532683*(RAL**5))+(0.20674845*(RAL**6)) 

+-(0.13264434*(RAL*V))-(0.00193383*(RAL**8)) 

ENDIF 

C 

IF  ((RAL  .GE.  0.55851)  .AND.  (RAL  .LT.  0.61087))  THEN 
C 

CYDAD=0.00023894+(0.00195121*(RAL-0.55851001))+(0.02459273 

+*(RAL-0.55851001)**2)-(0.1202244’((RAL-0.55851001)**3)) 

ENDIF 

C 

IF  (RAL  .GE.  0.61087)  THEN 
C 

CYDAD=0.27681 285-(2.02305395*RAL)+(6.01 1 8071 5*RAL*RAL) 

+-(9.242921 88*(R  AL**3))+(7.5985781 9*(RAL**4)) 
+-(2.8565527*(RAL**5))+(0.25460503*(RAL**7)) 

+-(0.01819815*(RAL*‘9)) 

ENDIF 

C 

C 

C  IF  (RAL  .LE.  0.0)  THEN 
C  EPA43=1.0d0 
C  ENDIF 

C  IF  (RAL  .GT.  0.0  AND  .LE.  0.6283185)  THEN 
’C  0.6283185  RADIANS  =  36  DEGREES 

C  EPA43=0.9584809+(4.13369452*RAL)-(18.31288396*RAL*RAL)+ 

C  +(19.551 1466*(RAL**3))-(1.09295946'RAL*DSPBD)+(0.17441033* 

C  +DSPBD*DSPBD) 

C  ENDIF 

C  IF  (RAL  .GT.,  0.6283185)  THEN 
EPA43=1.0cl0 
ENDIF 

•***«**#**•*•*«*****«*••«»•••****•»*•••*••»********•*«••»* 

*  NOTE  -  THE  PARAMETER  EPA43  IS  A  MULTIPLIER  ON  RUDDER 

*  EFFECTIVENESS  DUE  TO  SPEEDBRAKE.  THIS  TABLE  IS  ALSO 

*  LIMITED  TO  36  DEG  AOA.  HOWEVER,  THERE  IS  NO  AERODY 

*  NAMIC  EFFECT  FOR  ANGLES  OF  ATTACK  LESS  THAN  16  DEG, 

*  AND  THE  SPEEDBRAKE  IS  AUTOMATICALLY  RETRACTED  AT  AOA 

*  GREATER  THAN  15  DEG.  THEREFORE,  THIS  TABLE  SHOULD 
'  NOT  BE  NECESSARY  FOR  THE  ORDINARY  OPERATION  OF  THE 
'  AIRCRAFT 


CYDRD=0.00310199+(0.00119963'RAL)+(0.02806933*RAL*RAL) 
+-(0.1 2408447*(RAL**3))-(0.1 20321 21  *(RAL"4)) 
++(0.79150279*(RAL**5))-(0.86544347*(RAL**6)) 

++(0.278451 15*(RAL**7))+(0.00122999*RAL*RARUD)+(0.001 45943 
+*RARUD)-(0.01211427'RARUD*RARUD)+(0.00977937*(RARUD*‘3)) 

CYDTD=-0.001 57745-(0.0020881  *RAL)+(0.00557239'RAL'RAL) 
+-(0.00139886'(RAL**3))+(0.04956247*(RAL*'4)) 
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+-(0.0135353‘{RAL**5))-(0.11552397*(RAL**6)) 

++(0.11443452*(RAL**7))-(0.03072189*(RAL**8))-{0.01061113' 

+(RAL**3)*DELESRH0.00010529‘RAL*RAL'DELESR*DELESR) 

+-(0.00572463*RAL*DELESR‘DELESR) 

++(0.0l88536rRAL*RAL*DELESR)-(0.01412258*RAL*(DELESR*’3)) 

+-(0.00081776*DELESR)+(0.00404354*(DELESR*‘2))- 

+(0.00212189*(DELESR*‘3))+(0.00655063'(DELESR**4)) 

++(0.03341  584*(DELESR**5)) 

C 

RALY1  =0.61 08652 
RALY2=90.0d0/DEGRAD 
RBETY1  =-0.0872665 
RBETY2=0. 1745329 
C 

AY=0.1640d0 

ASTARY=0.95993 

BSTARY=0.087266 

C 

ZETAY={2.0D0'ASTARY-(RALY1 +RALY2))/(RALY2-RALY1 ) 
ETAY=(2.0D0‘BSTARY-(RBETY1+RBETY2))/(RBETY2-RBETY1) 

C 

X=(2.0D0*RAL-(RALY1+RALY2))/(RALY2-RALY1) 
Y=(2.0DO*RBETA-(RBETY1+RBETY2))/(RBETY2-RBETY1)  . 

C 

FY=({5.0DO*(2ETAY'‘2))-{4.0D0*ZETAY‘X)-1  .ODO)*(((X**2)-1  .ODO) 

+**2)*(1 .0D0/(((ZETAY**2)-1 .0D0)*’3)) 

C 

GY=((5.0D0'(ETAY"2))-(4.0D0'ETAY*Y)-1.0D0)*(((Y**2)-1.0D0)*'2) 

+*(1 .0D0/(((ETAY‘*2)-1 .0D0)*‘3)) 

C 

CYRB=AY*FY*GY 

C 

IF  (RAL  .LT,  0.6108652)  THEN 
C 

CYRB=O.OdO 
GOTO  500 
ENDIF 
C 

IF  ((RBETA  .LT.  -0.0872665)  .OR.  (RBETA  .GT.  0.1745329))  THEN 
C 

CYRB=0.0d0 
GOTO  500 
ENDIF 
C 

500  CFY=(CFY1  •EPA02L)+(CYDAD*DDA)+(CYDRD'DRUDD*DRFLX5*EPA43)+ 
+((CYDTD*DTFLX5)*DELEDD)+(CFYP'P8)+(CFYR'RB) 

++CYRB 
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DTFLX1=0.9750cl0 

DRFLX1=0.850dO 

C 

CML1=-0.00238235-(0.04616235*RAL)+(0.1 05531 68*RAL*RAL) 

++(0. 1 054 1 585*(RAL"3))-(0.40254765*(RAL**4)) 
++(0.3253049r(RAL**5))-(0.08496121*(RAL**6)) 

++(0.001 1 2288*(RAL**7))-(0.05940477*RABET*RAL)- 

+(0.07356236*RABET)-(0.00550119*RABET*RABET)+(0.00326191 

+*(RABET**3)) 

C 

IF  (RAL  .LT.  0.29671)  THEN 
C 

CMLP=-0.24963201-(0.03106297'RAL)+(0.1243063rRAL*RAL) 
+-(8.9527461 8*(RAL**3))+(1 00.331 09929*(RAL**4)) 
++(275.70069578*(RAL**5))-(1178.83425699'(RAL'*6)) 

+-(21 02.6681 1522*(RAL**7))+(2274.8978555r(RAL**8)) 

ENDIF 

C 

IF  ((RAL  .GE.  0.29671)  .AND.  (RAL  .LT.  0.34907))  THEN 
C 

CMLP=-0.1635261-(3.77847099*(RAL-0.29671001))+(147.47639465 
+'(RAL-0.29671 001  )**2)-(1 295.94799805*(RAL-0.29671 001  )**3) 
ENDIF 
C 

IF  (RAL  .GE.  0.34907)  THEN 
C 

CMLP=-1.37120291+(7.06112181*RAL)-(13.57010422*RAL*RAL) 

++(11.21323850*(RAL*'3)) 

+-(4.26789425'(RAL**4))+(0.623738r(RAL'*5)) 

ENDIF 

C 

IF  (RAL  .LT.  0.7854)  THEN 
C 

CMLR=0.03515391  +(0.59296381  *RAL)+(2.27456302*RAL*RAL) 
+-(3.8097803*(RAL**3)) 

+-(45.831  62842*(RAL'*4))+(55.31669213*(RAL-'5))+ 

+(194.29237485*(RAL**6))-(393.22969953*(RAL"7))+(192.20860739' 

+(RAL*'8)) 

ENDIF 

C 

IF  ((RAL  -GE.  0.7854)  .AND.  (RAL  .LE.  0.87266))  THEN 
C 

CMLR=0.0925579071-(0.6000000238'(RAL-0.7853999734)) 

++(1 .351 593971 3*((RAL-0.7a53999734)'*2)) 
++(29.0733299255*((RAL-0.7353999734)*'3)) 

ENDIF 

C 
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IF  (RAL  .GT.  0.87266)  THEN 
C 

CMLR=-311.126041+(1457.23391042*RAL)-(2680.19461944*RAL*RAL)+ 
+(2361 .4491 4738'(RAL**3))-{893.83567263*(RAL"4))+(68.23501 924* 
+(RAL**6))-(1.72572994*(RAL**9)) 

ENDIF 

G 

CLDAD=0.00057626+(0.00038479*RAL)-(0.00502091*RAL*RAL) 

++(0.001 61 407*(RAL**3))+{0.02268829*(RAL**4)) 
+-(0.03935269*(RAL**5))+(0.02472827‘(RAL**6)) 
+-(0.00543345*(RAL**7))+(0.0000007520348*DELESR*RAL)+ 
+(0.000000390773*DELESR) 

C 

CLDRD=0.00013713-(0.00035439*RAL)-(0.00227912*RAL‘RAL) 

++(0.00742636*(RAL**3))+(0.00991839*(RAL**4)) 

+-(0.0471 1846*(RAL**5))+(0.046124*(RAL**6)) 
+-(0.01379021*(RAL**7))+(0.00003678685*RARUD*RAL)+ 

+(0.00001 04375rRARUD)-(0.00015666*RARUD*RARUD)+(0.0001 61 33 
+*(RARUD**3)) 

G 

GLDTD=0.00066663+(0.00074174*RAL)+(0.00285735*RAL*RAL) 
+-(0.02030692*(RAL**3))-(0.00352997*(RAL**4)) 
++(0.0997962‘(RAL**5))-(0.1 4591 227* 

+(RAL**6))+(0.08282004*(RAL**7)) 

+-(0.01 68667*(RAL**8))+(0.003061 42*(RAL**3)*DELESR) 

+-(0.001 10266*RAL*RAL*(DELESR**2))+(0.00088031*RAL* 
+(DELESR**2))-(0.00432594*RAL*RAL*DELESR)- 
+(0.007201 41  *RAL*(DELESR**3)) 

+-(0.00034325*DELESR)+(0.00033433*(DELESR**2))+(0.00800183 

+*(DELESR**3))-(0.00555986*(DELESR**4))-(0.01841172*(DELESR**5)) 

G 

IF  (RAL  .LT.  0.0)  THEN 
G 

DGLB=-0.000060d0 

ENDIF 

G 

IF  ((RAL  .GE.  0.0)  .AND.  (RAL  .LE.  0.209434))  THEN 
G 

DGLB=-O.000060ci0+(0.0041035078*RAL*RAL)-(0.0130618699*(RAL**3)) 

ENDIF 

G 

IF  (RAL  .GT.  0.209434)  THEN 
G 

DGLB=0.0d0 

ENDIF 


GML=(GML1  *EPA02S)+(GLDAD*DDA)+(CLDRD*DRUDD*DRFLX1  •EPA43)+ 
+((GLDTD*DTFU1)*DELEDD)+(GMLP*PB)+(GMLR*RB)+(DGLB*BETA) 

G 
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*****************************  ******************1H**#************»»*1»0 


CMM1  =0.00501 496-{0.0800490rRAL)-(1.03486675*RAL*RAL) 
+-(0.68580677*(RAL**3))+(6.46858488*(RAL**4)) 
+-(10.155741uS*(RAL**5))+ 

+(6.44350808*(RAL*‘6))-(1 .461 751 88*(RAL**7)) 
++(0.24050902*RAL*DELESR) 

+-(0.42629958*DELESR)-(0.03337449*DELESR‘DELESR) 

+-(0.53951 733*(DELESR**3)) 

modified  25  Jul  90  to  use  new  curve  fit  for  CMQ 

OLD  EQUATION 

IF  (RAL  .LE.  0.25307)  THEN 

CMMQ=-3.8386262+(13.54661297*RAL)+(402.53011559*RAL*RAL) 

+-(6660.95327122*(RAL'*3))-(62257.89908743*(RAL*M') 

++(261 526. 1 0242329’(RAL**5)) 

++(2177190.33155227*(RAL**6))-(703575.13709062*(RAL**7))- 

+(20725000.34643054*(RAL*‘8))-(27829700.53333649*(RAL**9)) 

ENDiF 

IF  ((RAL  .GT.  0.25307)  .AND.  (RAL  .LT.  0.29671))  THEN 

CMMQ=-8.4926528931  -(2705.3000488281  '(RAL-0.2530699968)) 

++(1 23801 .5*(RAL-0.2530699968)**2) 
+-(1414377*(RAL-0.2530699968)**3) 

ENDIF 

IF  (RAL  .GE.  .29671)  THEN 

CMMQ=47.24676075-(709.60757056*RAL)+(3359.08807193’RAL*RAL)- 
+(7565.3201 7266*(RAL'*3))+(8695.1 858091  *(RAL*‘4)) 

+-(4891 .771 8331 3*(RAL“5))+(1 061 .5591 5089*(RAL"6)) 

ENDIF 

CMMQ  vs.  alpha  n  degrees 


NEW  EQUATION 

convert  alpha  to  degrees 

A=RAL*DEGRAD 
C 

FI  =-4.33509d0+A*(-0.1 41 624d0+A*(0.0946448d0+A*(-0.00798481  do 
+  +A*  (-0.001 68344dO+A*(0.000260037dO+A*(6.64054d-6+A*( 
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+  -2.20055d-6+A*(-2.74413d-8+A*(7.14476d-9+A* 

+  2.07046d-10))))))))) 
c 

F2=-302.567+a’(1 06.288+a*(-1 4.7034+A‘(1 .02524+ A*(-0.0393491 
+  +A*(0.00084082+A*(-9.365e-6+A*4.2355e-8)))))) 

C 

F3=1 724.99+ A*(-1 58.944+A*(5.59729+ A*(-0.0949624+ A*( 

+  0.000779066+A*(-2.47982e-6))))) 

C 

c  ramp  functions 
c 

R1=1.0-0.75*(A-10.0)**2+0.25*(A-10.0)**3 

R2=1.0-R1 

R3=1.0-7.5*(A-40.0)**2/62.5+(A-40.0)**3/62.5 

R4=1.0-R3 


C 


C 


IF(A.LT.10.0)THEN 

CMMQ=F1 

ELSEIF(A.LT.12.0)THEN 

CMMQ=FrR1+F2*R2 

ELSEIF(A.LT.40.0)THEN 

CMMQ=F2 

ELSEIF(A.LT.45.0)THEN 

CMMQ=F2*R3+F3*R4 

ELSE 

CMMQ=F3 

ENDIF 

CMM=CMM1+(CMMQ*QB) 


'*«*****«***«*«*«****•*****•***•*•***•«***«*****•** 


•**«»*******i»*i»* 


C 


DTFLX3=0.9750d0 

DRFLX3=0.890d0 


CMN 1  =0.01 441 51 2+(0.02242944*RALH0.30472558*(RAL**2)) 
++(0.14475549*(RAL**3)) 

++(0.931401 12*(RAL**4))-(1 .52168677*(RAL*'5))+ 

+(0.9074341 3*(RAL**6))-(0.1 651 0989*(RAL**7)) 

+-(0.0461 968*(RAL**8)) 

++(0.01 754292*(RAL**9))-(0.1 7553807‘RAL*RABET)+ 

+(0.15415649*RAL*RABET*DELESR) 

++(0.14829547*(RAL**2)*(RABET**2)) 

+-(0.1 1605031  ‘(RAL**2)*RABET*DELESR) 
+-(0.06290678*(RAL’'*2)'(DELESR**2)) 
+-(0.01404857*(RAL"2)*(DELESR**2)) 
++(0.07225609*RABET)-(0.08567087*(RABET**2)) 
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++(0.01 1 84674*(RABET**3)) 

+-(0.00519152*RAL*DELESR)+(0.03865177*RABET*DELESR) 

++{0.0006291 8*DELESR) 

C 

CNDRD=-0.001 53402+^0.001 84982*RAL)-{0.0068693*RAL*RAL) 
++(0.01772037*(RAL**3)) 

++(0.03263787*(RAL**4))-(0.15157163*(RAL**5))+(0.1 8562888 

+*(RAL**6))-(0.0966163*(RAL**7))+(0.01859168*(RAL**8))+(0.0002587 

+*RAL*DELESR)-(0.00018546*RAL*DELESR*RBETA)-(0.00000517304*RBETA) 

+-(0.00102718*RAL*RBETA)-(0.0000689379*RBETA*DELESR)-(0.00040536 

+*RBETA*RARUD)-(0.00000480484*DELESR*RARUD) 

+-(0.00041 786*RAL*RARUD) 

++(0.0000461872*RBETA)+(0.00434094*(RBETA**2)) 

+-(0.00490777*(RBETA**3)) 

++(0.0000051 57867*RARUD)+(0.00225169*RARUD*RARUD)-(0.00208072 
+*(RARUD*‘3)) 

C 

IF  (RAL  .LT.  0.55851)  THEN 
C 

CMNP=-0.00635409-(1 .141 53932*RAL)+(2.821 19027*(RAL**2))+ 
+(54.4739579*(RAL*’3))-(140.89527667*{RAL‘*4))-(676.73746128' 
+(RAL**5))+(2059.18263976*(RAL"6))+(1579.41664748*(RAL**7)) 
+-(8933.0853571 2*(RAL**8))+(6806.54761267*(RAL**9)) 

ENDIF 

C 

IF  ((RAL  .GE.  0.55851001)  .AND.  (RAL  .LE.  0.61087))  THEN 
C 

CMNP=-.07023239+(1 .08581 5*(RAL  -0.55851)) 
++(8.85265r((RAL-.55851)**2))-(192.6093‘((RAL-0.55851)“3)) 

ENDIF 

C 

IF  (RAL  .GT.  0.61087)  THEN 
C 

CMNP=-71 .03693533+(491 .3250671 5*RAL) 

+-(1388.1 1 1 77979’(RAL*'2))+ 

+(2033.48621905*(RAL**3)) 

+-(1 590.91 322362*(RAL**4))+(567.38432316*(RAL‘*5)) 
+-(44.97702536*(RAL'*7))+(2.8140669*(RAL**9)) 

ENDIF 


IF  (RAL  .LE.  -.069813)  THEN 
C 

CMNR=  -0.28050d0 
ENDIF 
C 

IF  ((RAL  .GT.  -.069813)  .AND.  (RAL  .LT.  0.0))  THEN 
C 

CMNR=-0.2804999948+(35.9903717041*(RAL+.0698129982)**2) 
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+-(51 6.1 574707031  *(RAL+ .06981 29982)**3) 

ENDIF 

C 

IF  ((RAL  .GE.  0.0)  .AND.  (RAL  .LE.  0.78539801))  THEN 
C 

CMNR=-.28071 511 -(2.521 83924*RAL)+(68.90860031*(RAL**2)) 
+-(573.2310051 1*(RAL**3))+(2009.08725005*(RAL*M)) 

+-(3385. 1 5675307*(R  AL**5)) 

++(2730.49473149‘(RAL**6))-(848.12322034*(RAL*‘7)) 

ENDIF 

C 

IF  ((RAL  .GT.  0.78539801)  .AND.  (RAL  .LT.  0.95993102))  THEN 
C 

CMNR=-0.1096954+(0.52893072*(RAL-0.78539801))-(6.09109497*(RAL- 
+0.78539801  )**2)+(1 7.4783401 5*(RAL-0.78539801  )**3) 

ENDIF 

C 

IF  (RAL  .GE.  0.95993102)  THEN 
C 

CMNR=-0.110d0 

ENDIF 

C 

CNDTD=0.00058286+(0.0007341*RAL)-(0.00746113*RAL*RAL) 

+-(0.00685223*(RAL**3)) 

++(0.03277271  *(RAL**4))-(0.02791456*(RAL**5)) 

++(0.0073291 5*(RAL**6)) 

++(0.001 20456*RAL*DELESR)-(0.00168102*DELESR)+(0.0006462* 
+DELESR‘DELESR) 

C 

CNDAD=0.00008228887-(0.00014015*RAL)-(0.0013493*RAL*RAL)+ 

+(0.00020487*(RAL**3))+(0.0056124r(RAL**4)) 

+-(0.00634392*(RAL**5)) 

++(0.00193323*(RAL**6))-(2.05815E-17*(RAL*DAILA))+(3.794816E-17' 

+(DAILA‘*3)) 

C 

DCNB=-2.500E-4 

C 

RALN1  =0.6981 3 
RALN2=90.bd0/DEGRAD 
RBETN1  =-0.1 74532 
RBETN2=0.34906 
C 

AN=0.034d0 

ASTARN=1.0472d0 

BSTARN=0.087266 

C 

ZETAN=(2.0D0'ASTARN-(RALN1+RALN2))/(RALN2-RALN1) 

ETAN=(2.0D0*BSTARN-(RBETN1+RBETN2))/(RBETN2-RBETN1) 

C 

XN=(2.0D0*RAL-(RALN1+RALN2))/(RALN2'RALN1) 
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YN=(2.0D0*RBETA-(RBETN1+RBETN2))/(RBETN2-RBETN1) 

C 

FN=((5.0D0*{ZETAN‘’2))-(4.0D0*ZETAN*XN)-1.0D0)* 
+(((XN**2)-1 .0D0)'*2)/(((ZETAN**2)-1 .0D0)**3) 

C 

GN=((5.OD0*(ETAN**2))-(4.OD0*ETAN*YN)-1.OD0)* 

+(((YN**2)-1 .0D0)**2)/(((ETAN**2)-1 .0D0)**3) 

C 

CNRB=AN*FN*GN 

C 

IF  (RAL  .LT.  0.69813)  THEN 
C 

CNRB=0.0d0 

GOTO  1000 

ENDIF 

C 

IF  ((RBETA  .LT.  -0.174532)  .OR.  (RBETA  .GT.  0.34906))  THEN 
C 

CNRB=0.0d0 

GOTO  1000 

ENDIF 


1000  CMN=(CMNrEPA02S)+(CNDAD'DDA)+((CNDRD*DRUDD*DRFLX3)*EPA43)+ 
.  +((CNDTD*DTFLX3)*DELEDD)+(CMNP*PB)+{CMNR*RB)+(DCNB*BETA) 
++CNRB 
C 
C 

^*«  ***************  ****•'«*********#*******•  ******«****«***********<»****^ 

C 

C  THIS  SECTION  DETERMINES  THE  EFFECT  OF  THE  THRUST  VALUES  FOR 
C  ADDITION  TO  CX,  CY,  CZ,  CLM,  CMM,  AND  CNM  VALUES  DETERMINED 
C  ABOVE  AND  CONTAIN  THE  FOLLOWING  VARIABLES: 

C  CPTAL  -  COSINE  OF  PITCH  VECTOR  ANGLE 
C  SPTAL  -  SINE  OF  PITCH  VECTOR  ANGLE 
C  CYTAL  -  COSINE  OF  YAW  VECTOR  ANGLE 
C  SYTAL  -  SINE  OF  YAW  VECTOR  ANGLE 
C  ENGPQ  -  PORT  ENGINE  THRUST/(QBAR*S) 

C  ENGSQ  -  STARBOARD  ENGINE  THRUST/(QBARS*S) 

C  CXENGP  -  COEFFICIENT  OF  PORT  ENGINE  THRUST  IN  X  DIRECTION 

C  CXENGS  -  COEFFICIENT  OF  SBRD  ENGINE  THRUST  IN  X  DIRECTION 

C  CXT  -  COEFFICIENT  OF  TOTAL  THRUST  IN  X  DIRECTION 
C  CYENGP  -  COEFFICIENT  OF  PORT  ENGINE  THRUST  IN  Y  DIRECTION 

C  CYENGS  -  COEFFICIENT  OF  SRBD  ENGINE  THRUST  IN  Y  DIRECTION 

C  CYT  -  COEFFICIENT  OF  TOTAL  THRUST  IN  Y  DIRECTION 
C  CZENGP - 
C  CZENGS - 
C  CZT  - 
C  CLMT  - 
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C  CMMT  - 
C  CNMT  - 
C 

CPTAL=COS(PTAL) 

SPTAL=SIN(PTAL) 

CYTAL=COS(YTAL) 

SYTAL=SIN(YTAL) 

CRAL=COS(RAL) 

SRAL=S1N(RAL) 

C 

ENGPQ=ENGP/QBARS 

ENGSQ=ENGS/QBARS 

C 

CXENGP=ENGPQ*CPTAL*CYTAL 

CXENGS=ENGSQ'CPTAL*CYTAL 

CXT=CXENGP+CXENGS 

C 

CYENGP=ENGPQ'CPTAL‘SYTAL 

CYENGS=ENGSQ'ePTAL*SYTAL 

CYT=CYENGP+CYENGS 

C 

C2ENGP=ENGPQ*SPTAL 

CZENGS=ENGSQ*SPTAL 

CZT=CZENGS+C2ENGP 

C 

CLMT=(CZENGS-eZENGP)*(25.6d0/12.0d0)/BWING 

C 

CMMT=CXT*(0.25d0/1 2.0d0)/CWING+ 

+  CZT*20.219d0/CWING 
C 

CNMT=:(CXENGP-CXENGS)'(25.5d0/12.0d0)/BWING- 
+  CYT*20.219d0/BWING 
C 

CX=CFZ‘SRAL-CFX*CRAL+CXT 

CY=CFY+CYT 

CZ=-(CFZ*CRAL+CFX*SRAL)+CZT 

CLM=CML+CLMT 

CMM=CMM+CMMT 

CNM=CMN+CNMT 


THE  0.25/12.0  IS  THE  OFFSET  OF  THE  THRUST  VECTOR  FROM  THE  CG 
THE  20.219  is  the  moment  arm  from  the  nozzle  pivot  to  the  eg 
THE  25.5/12.0  is  the  moment  arm  of  the  engines  from  the  eg 

RETURN  CX,  CY,  CZ,  CLM,  CMM,  CNM  TO  CALLING  PROGRAM. 


RETURN 

END 
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THE  FOLLOWING  SECTIONS  ADDED  13  AUG  91 


********************************************************* 

THIS  SUBROUTINE  EXECUTES  THREE  DIFFERENT  BLENDS  FROM 
THE  MCDONNELL  MODEL  COEFICIENT  DATABASE  TO  THE  ROTARY 
BALANCE  COEFICIENTS  GIVEN  THE  VALUE  OF  PARAMETER  (10); 

BLEND.  THE  FIRST  MODE  (PAR(10)<1.0)  IS  INTENDED  TO 
BE  USED  TO  DETERMINE  AN  EQUILIBRIUM  STATE  SET  OF 
PARAMETERS  FOR  THE  ROTARY  BALANCE  DATA  MODEL  GIVEN  A 
SET  OF  KNOWN  STATIC  AERODYNAMIC  EQUILIBRIA.  THE 
SECQND  MODE  (PAR(IO)  =  2.0)  SETS  THE  MODEL  TO  EXECUTE 
PURE  ROTARY  BALANCE  DATA  ABOVE  AOA  35  DEGREES.  THE 
THIRD  MODE  (PAR(IO)  =  3.0)  IS  A  HYBRID  MODEL  WHERE 
BOTH  THE  STATIC  AND  ROTARY  BLANCE  DATABASES  ARE 

COMBINED  UTILIZING  UNIQUE  ASPECTS  OF  BOTH  DATABASES. 
>****«************«************<************************** 


SUBROUTINE  RBBLEND(U,PAR,NDIM.ICP) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  /SEIZE/  CX.CY.CZ.CLM.CMM.CNM 
COMMON  /SEIZET/  CXT,CYT,CZT,CLMT,CMMT,CNMT 
COMMON  /SEIZER/  CXR,CYR,CZR,CLMR,CMMR,CNMR 
COMMON  /RBPOLY/  PC,  OMEGA 
DIMENSION  PAR(1 0),PC(6.79),U(NDIM) 


IF  (PAR(IO)  .LT.  2.0)  THEN 


BLEND  =  PAR(IO) 

C 

CXS  =  CX-CXT 
CYS  =  CY-CYT 
CZS  =  CZ-CZT 
CLMS  =  CLM-CLMT 
CMMS  =  CMM-CMMT 
CNMS  =  CNM-CNMT 
C 

CX  =  CX  +  (CXR-CXS)*BLEND 
CY  =  CY  +  (CYR-CYS)'BLEND 
CZ  =  CZ  +  (CZR-CZS)'BLEND 
CLM  =  CLM  +  (CLMR-CLMS)*BLEND 
CMM  =  CMM  +  (CMMR-CMMS)*BLEND 
CNM  =  CNM  +  (CNMR-CNMS)*BLEND 
C 

ENDIF 
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IF  {PAR(IO)  .EQ.  3.0)  THEN 
C 

AL=  U(1) 

BETA  =  U{2) 

DEGRAD  =  57.2957795 
RAL  =  AL/DEGRAD 
RBETA  =  BETA/DEGRAD 
C 

PTEMP  =  U(3) 

QTEMP  =  U(4) 

RTEMP  =  U(5) 

C 

P  =  U(3)  -  OMEGA*DSIN(RBETA) 

Q  =  U(4)  -  OMEGA*'DSIN(RAL)*DCOS(RBETA) 
R  =  U(5)  -  OMEGA*DCOS(RAL)*DCOS(RBETA) 
C 

U(3)  =  P 
U(4)  =  Q 
U(5)  =  R 
C 

CXRP  =  CXR 
CYRP  =  GYR 
CZRP  =  CZR 
CLMRP  =  CLMR 
CMMRP  =  CMMR 
CNMRP  =  CNMR 
C 

CALL  COEFF(U,PAR,NDIM,ICP) 

C 

U(3)  =  0 
U{4)  =  0 
U(5)  =  0 
C 

CALL  RBCOEF(U,PAR,NDIM) 

C 

CXRO  =  CXR 
CYRO  =  CYR 
CZRO  =  CZR 
CLMRO  =  CLMR 
CMMRO  =  CMMR 
CNMRO  =  CNMR 
C 

CXS  =  CX-CXT 
CYS  =  CY-CYT 
CZS  =  CZ-CZT 
CLMS  =  CLM-CLMT 
CMMS  =  CMM-CMMT 
CNMS  =  CNM-CNMT 
C 
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CXR  =  CXS  +  (CXRP-CXRO) 

CYR  =  CYS  +  (CYRP-CYRO) 

CZR  =  CZS  +  (CZRP-CZRO) 

CLMR  =  CLMS  +  (CLMRP-CLMRO) 
CMMR  =  CMMS  +  (CMMRP-CMMRO) 
CNMR  =  CNMS  +  (GNMRP-CNMRO) 
C 

U{3)  =  PTEMP 
U(4)  =  QTEMP 
U(5)  =  RTEMP 
C 

CALL  COEFF(U,PAR,NDIM,ICP) 

C 

ENDIF 


IF  (PAR(IO)  .GE.  2.0)  THEN 


DEGRAD  =  57.2957795 
AL  =  U(1) 

RAL  =  AL/DEGRAD 

DELTA  =  5/DEGRAD 

BLEND  =  (RAL  -  0.523598776)/DELTA 

IF  ((RAL  .LE.  0.610865238)  .AND.  (RAL  .GE.  0.523598776))THEN 

CX  =  CX  +  (CXR-CX+CXT)*(3-2*BLEND)*BLEND*BLEND 

CY  =  CY  +  (CYR-CY+CYT)*(3-2*BLEND)*BLEND*BLEND 

CZ  =  CZ  +  (CZR-CZ+CZT)*(3-2*BLEND)*BLEND*BLEND 

CLM  =  CLM  +  (CLMR-CLM+CLMT)*(3-2*BLEND)*BLEND‘BLEND 

CMM  =  CMM  +  (CMMR-CMM+CMMT)*(3-2*BLEND)*BLEND-BLEND 

CNM  =  CNM  +  (CNMR-CNM+CNMT)'(3-2'BLEND)*BLEND’BLEND 

ELSE 

CX  =  CXR  +  CXT 
CY  =  CYR  +  CYT 
CZ  =  CZR  +  CZT 
CLM  =  CLMR  +  CLMT 
CMM  =  CMMR  +  CMMT 
CNM  =  CNMR  +  CNMT 

ENDIF 

ENDIF 

RETURN 

END 
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*  THE  FOLLOWING  TV/O  SUBROUTINES,  R8POLYCOEF  AND  RBCOEF, 

*  REPRESENT  THE  ROTARY  BALANCE  COEFICIENT  DATA.  THE  FIRST 

*  ROUTINE  READS  IN  THE  COEFIGIENTS  OF  THE  POLYNOMIALS  THAT 

*  REPRESENT  THE  SIX  COEFIGIENTS  AND  THE  SECOND  ROUTINE  IS  THE 

*  POLYNOMIAL  THAT  REPRESENTS  THE  POLYNOMIAL 


********************** ********** ************ ***ftft*****************i^ 

C 

SUBROUTINE  RBPOLYCOEF 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-2) 

COMMON  /RBPOLY/  PC,  OMEGA 
DIMENSION  PC(6,79) 

C 

0PEN(99,FILE=’RBDATA.150’,STATUS=’0LD’) 

REWIND(99) 

DO  110  1=1,6 
DO  100  L=1,79 
READ(99,*)  PC(I,L) 

100  CONTINUE 
110  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  RBCOEF(U,PAR,NDIM) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON  /SEI2ER/  CXR,CYR,CZR,CLMR,CMMR,CNMR 
COMMON  /RBPOLY/  PC,  OMEGA 
DIMENSION  U(NDIM),PAR(1 0),PC(6,79) 

C 

AL  =U(1) 

DEGRAD  =  57.29577951 
C 

RAL  =  AL/DEGRAD 
C 

IF  ((RAL  .LE.  1.5708)  .AND.  (RAL  .GE.  0.5235988))  THEN 
C 

CALL  RBDATA(U,PAR,NDIM,1,CXR) 

CALL  RBDATA(U,PAR.NDIM,2,CYR) 

CALL  RBDATA(U,PAR,NDIM,3,CZR) 

CALL  RBDATA(U,PAR,NDIM,4,CLMR) 

CALL  RBDATA(U,PAR,NDIM,5,CNMR) 

CALL  RBDATA(U,PAR,NDIM,6,CMMR) 

CXR  =  -CXR 
CZR  =  -CZR 
ENDIF 
RETURN 
END 
C 
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SUBROUTINE  RBDATA(U, PAR, NDIM, I, CF) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
COMMON/ACDATA/BWING,CWING.SREF.RHO,RMASS 
COMMON  /RBPOLY/  PC,  OMEGA 
DIMENSION  U(NDIM),PAR(1 0),PC(6,79) 

C 

DEGRAD  =  57.29577951 
C 

A  =  U(1)/DEGRAD 
B  =  U(2)/DEGRAD 
P  =  U(3) 

0  =  U(4) 

R  =  U(5) 

C 

VTRFPS  =  U(8)*1000 
C 

DE  =  PAR(1)/DEGRAD 
DR  =  PAR(2)/DEGRAD 
DA  =  PAR(3)/DEGRAD 
DD  =  PAR(9)/DEGRAD 

FOR  COMPARISON  TO  PREVIOUS  MODELS,  DIFFERENTIAL  ELEVATOR 
WILL  BE  VARIED  AS  A  FUNCTION  OF  AILERON  DEFLECTION. 

DD  =  0.3*DA 

OMEGA  =  (P*OCOS(A)  +  R‘DSIN(A))*DCOS{B)  +  Q*DSIN(B) 

OM  =  (BWING*OMEGA)/(2'VTRFPS) 

79  TERM  ROTARY  BALANCE  DATA  POLYNOMIAL 


Cl  =  PC(I,1)  +  PC(I.2)*A  +  PC(I,3)'B 

C2  =  PC(I,4)'DA  +  PC(I.5)'DD  +  PC{I,6)*DE 

C3  =  PC(i,7)*OM  +  PC(I,8)'DR  +  FC(I,9)'A*A 

C4  =  PC(I,10)'DA*B  +  PC(I,11)'B*B  +  PC(I.12)'B*DR 

C5  =  PC(l,13)'A'OM  +  PC(I,14)'DA*DD  +  PC(I,15)*DR*DR 

C6  =  PC(l,16)*OM*OM  +  PC(I,17)'DA*A  +  PC(I,18)*B*0M 

C7  =  PC(I,19)'DD*DD  +  PC(I,20)*OM-DE  +  PC(I,21)*DD*A 

C8  =  PC{I,22)'A*DE  +  PC(I,23)'B*DE 

C9  =  PC{I,24)*B*A  +  PC(l,25)*OM*DD  +  PC(I.26)*B*DD 

CIO  =  PC(l,27)*OM'DR  +  PC(l.28)*DA*OM'OM  +  PC(I,29)*DA'A‘A 

C11  =  PC(I,30)*DD*DD*B+PC(I,31)'OM*OM*A+PC(I,32)'OM'OM*OM 

C12  =  PC{I,33)'A'B*B  +  PC(l,34)*B'A'OM 

C13  =  PC{I,35)*B*A*DR  +  PC(I,36)*A-B'DD  +  PC(l,37)'B*B‘OM 

C14  =  PC(l,38)*DD*OM*OM  +  PC(l,39)'B'DA*OM 

C15  =  PC(I,40)*B'OM*DR 

C16  =  PC(I,41)*B*DA*A  +  PC(I,42)*DD*DD*A 

C17  =  PC(I,43)*A*DR'DR 

C18  =  PC(I,44)*DD'A'A  +  PC(I,45)*A*A-DE  +  PC(I,46)-B*A*A 
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C19  =  PC(l,47)*B*OM*DE+PC(l,48)*OM‘DR*DR+PC(l,49)*B‘OM*OM 
C20  =  PC(I,50)*A*A*OM 
C21  =  PC(1.51)*OM*OM*DE+PC(l,52)*A*OM‘DR 
C22  =  PC(l,53)*DD*DD*OM  +  PC(l.54)*A*A*A+PC{l,55)*A*DD*OM 
C23  =  PC(l.56)*A*A*DR+PC(l,57)*OM*OM*DR+PC(l,58)*A*OM*DE 
C24  =  PC(l,59)*A*A*OM*OM  +  PC(I.60)*B*B‘OM*OM 
C25  =  PC(I.61)*DD*DD*0M*0N/1 
C26  =  PC(l.62)*bD*DD*DD*A  +  PC(I,63)*B*DD‘DD*DD 
C27  =  PC(I.64)*A'A*B*B  +  PC(l,65)*B*OM*OM*OM 
C28  =  PC(l,66)*OM*OM*OM*OM  +  PC(l,67)*OM*OM*OM*DE 
C29  =  PC(l,68)*A*OM*OM*OM  +  PC(I,69)*DD*DD*A*A 
C30  =  PC(I,70)*DE*A*A*A  +  PC(I,71)*DA*A*A*A 
C31  =  PC(I,72)*B*A*A*A  +  PC(I,73)*A‘A*A*A 
C32  =  PC(I.74)*A'A'A*DD  +  PC(l,75)*OM*DD*DD*DD 
C33  =  PC(I,76)*A*A*A*DR  +  PC(l,77)*A‘A*DR*bR 
C34  =  PC(l,78)*OM*OM*DR*DR  +  PC(l,79)»OM*OM*OM*DR 
C 

CFl  =  CI+C2+  C3  +  C4  -T  C5  +  C6  +  C7  +  C8  +  C9  +  CIO 

CF2  =  Gil  +  C12  +  C13  +  C14  +  C15  +  C16  +  Cl7  +  CI8 

CF3  =  C19  +  C20  +  C21  +  C22  +  C23  +  C24  +  C25  +  C26 

CF4  =  C27  +  C28  +  C29  +  C30  +  C31  +  C32  +  C33  +  C34 

C 

CF  =  GF1  +  CF2  +  GF3  +  CF4 
C 

RETURN 

END 
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